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ABSTRACT 

-This  report  describes  in  detail  a  computational  method  for 
the  solution  of  the  partially-parabolic  (or  semi-elliptic  or 
parabolized)  Reynolds-averaged  Navier-Stokes  equations  for  the 
calculation  of  external  flow  around  ship-like  three-dimensional 
bodies.  Among  the  main  features  of  the  method  are  the  following. 
Numerically-generated  body-fitted  coordinates  are  used  to  faci¬ 
litate  applications  to  a  wide  variety  of  shapes.  The  convective- 
transport  equations  are  discretized  using  the  finite-analytic 
scheme  which  employs  analytic  solutions  of  the  locally-linearized 
equations.  A  time-marching  algorithm  is  employed  to  enable 
future  extensions  to  be  made  to  handle  unsteady  and  fully- 
elliptic  problems.  A  two-step  global  pressure-correction  algor¬ 
ithm  has  been  developed  to  accelerate  convergence.  The  method 
can  be  used  with  large  solution  domains  in  order  to  capture  the 
viscous-inviscid  interaction  so  that  iterative  matching  between 
separate  viscous-flow  and  potential-flow  solutions  is  not 
necessary.  For  turbulent  flows,  the  well  known  k^e  model  is 
used,  but  with  a  more  convenient  and  realistic  treatment  of  the 
flow  close  to  solid  walls.  < 

The  generality  and  flexibility  of  the  method  are  demon¬ 
strated  by  applications  to  several  two-dimensional,  axisymmetric 
and  three-dimensional  flows. 
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CALCULATION  OF  TRAILING-EDGE ,  STERN  AND  HAKE  FLOWS 
BY  A  TIME-MARCHING  SOLUTION  OF  THE  PART I ALLY—  PARABOL I C  EQUATIONS 


I.  INTRODUCTION 

Numerical  solutions  of  the  complete  Navier-Stokes  equations 
for  laminar  flow,  and  the  corresponding  Reynolds-averaged  equa¬ 
tions  for  turbulent  flow,  have  received  a  great  deal  of  attention 
in  recent  years  since,  in  principle,  they  describe  flows  with  any 
level  of  complexity,  the  only  uncertainty  being  that  introduced 
by  the  turbulence  model  employed  to  effect  closure  of  the  Rey¬ 
nolds  equations.  Indeed,  solutions  using  a  variety  of  finite- 
difference,  finite-element  and  hybrid  schemes  are  becoming  quite 
common  for  flows  involving  simple  two-dimensional  and  axisym- 
metric  geometries  or  some  other  simplifying  features.  For 
example,  many  more  solutions  have  been  obtained  for  internal 
flows  in  which  the  boundary  conditions  are  easier  to  specify. 
Due  to  the  explosion  of  the  literature  in  this  branch  of  computa¬ 
tional  fluid  mechanics  it  is  not  surprising  that  an  authoritative 
review  of  the  state-of-the-art  has  not  yet  appeared.  No  such 
review  is  intended  here.  However,  it  is  of  interest  to  note  here 
that  considerable  confusion  abounds  in  the  terminology  used  in 
the  literature,  particularly  with  regard  to  the  use  of  the  term 
"Navier-Stokes”.  In  addition  to  the  complete  Reynolds-averaged 
equations  for  turbulent  flow,  it  is  often  associated,  rather 
loosely  and  incorrectly,  with  equations  in  which  further  approxi¬ 
mations  are  made.  This  distinction,  or  the  lack  of  it,  is  par¬ 
ticularly  evident  in  solution  methods  which  claim  to  address 
three-dimensional  external  flows.  If  the  correct  definition  is 
insisted  upon,  then  there  exist  a  very  few  solutions  of  the 
Navier-Stokes  equations  or  the  Reynolds-averaged  equations  in 
three-dimensional  external  flows.  Methods  which  appear  to  have 
this  capability  are  those  of  Briley  and  McDonald  (1977), 


Spradley,  Stalnaker  and  Ratliff  (1981),  and  Hankey,  Graham  and 
Shang  (1982)  for  laminar  flow,  and  Shang  et  al.  (1980,  1983)  and 
Spalding  (1981)  for  turbulent  flow. 

As  the  power  of  modern  computers  continues  its  dramatic 
increase  and  the  costs  continue  to  decrease,  many  problems  of 
practical  interest,  which  invariably  involve  complex  three-dimen¬ 
sional  geometries,  will,  in  the  long  run,  become  tractable 
through  the  solutions  of  the  complete  equations.  The  rapid  pro¬ 
gress  being  made  in  another  branch  of  computational  fluid  mechan¬ 
ics,  namely  the  numerical  generation  of  computat ional  grids  for 
arbitrary  bodies,  will  undoubtedly  facilitate  such  applica¬ 
tions.  However,  the  ultimate  success  in  the  prediction  of  turbu¬ 
lent  flows  will  be  determined,  in  large  measure,  by  future  pro¬ 
gress  in  the  area  of  turbulence  modeling.  For  example,  compari¬ 
sons  made  at  the  1980-81  Stanford  Conferences  (Kline,  Cantwell 
and  Lilley,  1983)  for  simple  flows  involving  separation  indicated 
that  the  solutions  of  the  Reynolds-averaged  Nav ier-Stokes  equa¬ 
tions  show  considerable  dependence  on  the  turbulence  model 
used.  Ttius ,  the  modeling  of  turbulence  remains  the  major  source 
of  uncertainty.  Another  major  obstacle  to  progress  in  the  treat¬ 
ment  of  complex  three-dimensional  flows  is  the  difficulty  (and 
expense)  of  obtaining  comprehensive  and  reliable  data  which  can 
be  used  to  validate  the  solutions. 

The  primary  difficulty  in  the  solution  of  the  full  Navier- 
Stokes  equations  is  that  they  are  elliptic  and  consequently  the 
solution  must  be  found  simultaneously  in  all  three  spatial  direc¬ 
tions.  The  computer-storage  requirements  and  the  computational 
effort  can  be  greatly  reduced,  however,  by  the  use  of  equations 
whose  complexity  lies  somewhere  between  that  of  the  Navier-Stokes 
and  the  boundary-layer  equations.  Two  different  approaches  have 
been  followed  in  the  formulation  of  such  "higher-order"  equa¬ 
tions:  one  seeking  to  simplify  the  Navier-Stokes  equations  by 
discarding  certain  terms  and  the  other  to  generalize  the  first- 
order  boundary- layer  equations  by  introducing  additional  terms. 
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Although  the  two  approaches  lead  to  the  same  equations  in  the 
limit,  intermediate  equations  representing  different  levels  of 
approximations  have  been  used,  and  appropriate  solution  proce¬ 
dures  developed,  for  the  description  of  restricted  classes  of 
f lows . 

For  general  three-dimensional  flows,  we  can  distinguish  at 
least  three  levels  of  approximations  between  the  classical  boun¬ 
dary-layer  equations  and  the  complete  Navier-Stokes  equations.  A 
set  of  second-order  boundary-layer  equations  can  be  deduced  by 
formally  retaining  the  second-order  terms  in  the  familiar  order- 
of-magnitude  analysis  (Nash  and  Patel,  1972).  These  equations 
contain  additional  terms  involving  the  surface  curvatures  and  the 
variations  of  the  coordinate  metrics  in  the  direction  normal  to 
the  surface  which  are  neglected  in  the  first-order  equations. 
However,  the  approximations  inherent  in  such  an  analysis  restrict 
the  resulting  equations  to  shear  layers  developing  parallel  to 
continuously  curved  surfaces.  A  greater  degree  of  generality  is 
embodied  in  the  somewhat  inappropriately  named  'thin-layer'  or 
'thin-layer  Navier-Stokes'  equations,  in  which  all  viscous  and 
turbulent  transport  terms  other  than  those  in  the  direction  nor¬ 
mal  to  the  surface  are  neglected.  The  difference  between  these 
and  the  second-order  equations  lies  in  the  retention  of  the  con¬ 
vective  terms  in  the  normal  direction.  It  is  clear  that  both  are 
restricted  to  flows  in  which  cross-stream  diffusion  parallel  to 
the  surface  is  negligible  and  therefore  do  not  apply,  for  exam¬ 
ple,  to  the  flow  in  a  streamwise  corner.  The  next  level  of  gen¬ 
erality,  beyond  which  lie  the  complete  equations,  is  achieved  by 
the  so-called  part ially-parabol ic  (or  'parabolized  Navier- 
Stokes',  or  'semi-elliptic  Navier-Stokes')  equations,  in  which 
only  the  longitudinal  transport  due  to  viscosity  and  turbulence 
is  neglected. 

All  three  sets  of  equations  share  a  common  simplification, 
namely  the  destruction  of  the  ellipticity  of  the  Navier-Stokes 
equations  in  the  longitudinal  direction.  This  enables  the  solu- 
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tions  to  be  marched  in  the  downstream  direction,  from  some  ini¬ 
tial  plane,  provided  of  course  there  is  no  recirculation  or  re¬ 
verse  flow  in  that  direction.  Physically,  these  approximations 
differ  from  those  in  first-order  boundary-layer  theory  insofar  as 
the  pressure  is  regarded  as  an  unknown  and  provides  the  mechanism 
for  the  propagation  of  information  from  downstream.  If  the  down¬ 
stream  effects  are  weak,  as  in  thin  boundary  layers  and  in  fully- 
developed  flows  in  straight  ducts  of  arbitrary  cross-section, 
then  the  problem  is  truly  parabolic  and  a  single  downstream-mar¬ 
ching  solution  suffices.  If  they  are  not,  as  in  the  flow  over 
the  stern  and  in  the  near  wake  of  a  ship,  several  sweeps  are 
necessary  to  obtain  the  solution.  Methods  utilizing  the  more 
general  equations  differ  in  the  order  in  which  the  continuity  and 
momentum  equations  are  solved  and  in  the  manner  in  which  the 
pressure  is  updated  in  each  sweep. 

It  should  be  noted  that  the  thin-layer  and  partially-para- 
bolic  equations  are  identical  in  two-dimensional  and  axisymmetric 
flows  due  to  the  absence  of  transverse  diffusion.  The  substan¬ 
tial  simplification  introduced  by  two  dimensionality  has  led  to 
the  development  of  several  different  numerical  methods  for  the 
solution  of  these  equations.  Some  of  these  will  be  discussed 
later. 

In  three-dimensional  flows,  the  thin-layer  approximations 
were  introduced  by  Baldwin  and  Lomax  (1978)  to  treat  a  class  of 
flows  in  high-speed  aerodynamics.  The  most  recent  applications 
are  those  of  Pulliam  and  Steger  (1980),  Hung  (1980)  and  Chien  and 
Hsieh  (1983)  for  laminar  flow,  and  Pulliam  and  Lomax  (1979), 
Deiwert  (1981)  and  Hung  and  Chaussee  (1981)  for  turbulent  flow. 
These  approx ima t ions  were  introduced  primarily  to  reduce  the 
computational  effort  involved  in  the  solution  of  the  complete 
Nav ier-Stokes  equations  and,  as  pointed  out  by  Hung  and  Chaussee, 
the  thin-layer  solutions  for  supersonic  flow  past  bodies  at  inci¬ 
dence  agree  with  those  of  the  Nav ier-Stokes  equations.  Thin- 
layer  equations  have  also  been  used,  quite  independently,  in  ship 
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hydrodynamics.  However,  they  have  been  obtained  by  an  entirely 
different  route,  namely  by  the  generalization  of  the  boundary- 
layer  equations  to  account  tor  the  thick  boundary-layer  effects 
over  ship  sterns.  The  equations  used  by  Soejima  (1983)  and  Od- 
abasi  and  Davies  (1983),  and  the  corresponding  higher-order  mo¬ 
mentum-integral  equations  of  Nagamatsu  (1980)  and  Larsson  and 
Chang  (1980),  belong  in  this  latter  category.  Although  the  un¬ 
derlying  assumptions  are  similar,  and  lead  to  identical  equations 
for  incompressible  flow,  the  solution  methodology  has  also  deve¬ 
loped  along  different  lines  since  the  supersonic-flow  applica¬ 
tions  in  aerodynamics  are  viewed  from  a  Nav ier-Stokes  perspective 
whereas  those  in  ship  hydrodynamics  continue  to  be  viewed  from  a 
boundary-layer  perspective.  While  the  former  have  been  quite 
successful,  the  latter  have  met  with  only  limited  success,  pre¬ 
sumably  because  extensions  of  solution  procedures  developed  for 
thin  boundary  layers  may  not  be  the  most  appropriate  for  the  more 
general  equations.  A  major  difference  between  the  two  approaches 
lies  in  the  manner  in  which  the  important  upstream  influence  of 
the  pressure  field  is  treated. 

As  noted  earlier,  the  part ial ly-parabol ic  equations  differ 
from  the  thin-layer  equations  in  that  terms  in  viscous  and  turbu¬ 
lent  transport  of  momentum  in  the  transverse  direction  parallel 
to  the  surface  are  retained.  Apart  from  the  problem  of  modeling 
the  additional  Reynolds  stresses,  there  is  no  difference  in  the 
mathematical  properties  of  the  two  sets  of  equations,  and  there¬ 
fore,  essentially  the  same  numerical  techniques  can  be  used  for 
their  solution.  The  partially-parabolic  approximations  were 
first  introduced  by  Pratap  and  Spalding  (1976),  to  describe  flows 
in  which  there  is  a  predominant  flow  direction  (along  which  the 
diffusion  is  neglected  and  there  is  no  recirculation)  but  signi¬ 
ficant  pressure  effects  propagate  upstream.  The  retention  of  the 
transverse  diffusion  terms  enables  the  treatment  of  a  much  wider 
class  of  three-dimensional  flows  of  practical  interest,  e.g.  the 
flow  in  curved  ducts  of  arbitrary  cross-section,  the  flow  in  a 


streamwise  corner  (wing-fuselage  or  hull-appendage  junctions), 
the  flow  over  the  stern  and  in  the  wake  of  ship,  and  the  flow 
over  slender  bodies  at  incidence. 

The  solution  of  the  part ially-parabol ic  equations  in  three- 
dimensional  external  flows  is  the  subject  of  the  present 
report.  This  report  describes  the  development  of  a  numerical 
method  which  incorporates  several  novel  features  designed  to 
improve  the  accuracy  of  the  solutions  and  the  ease  of  application 
to  a  wide  range  of  flows  of  practical  interest.  The  types  of 
flows  which  can  be  treated  by  the  part ial ly-parabol ic  equations 
are  described  in  the  next  section  along  with  the  general  features 
of  the  new  method.  Subsequent  sections  are  concerned  with  the 
transformation  of  the  equations  into  curvilinear  body-fitted 
coordinate  systems,  the  numerical  generation  of  the  coordinates 
and  the  novel  finite-analytic  method  used  to  solve  the  equa¬ 
tions.  The  report  concludes  with  the  presentation  of  some  solu¬ 
tions  in  two-dimensional,  ax isymmetric,  and  three-dimensional 
flows  which  demonstrate  the  performance  and  capabilities  of  the 
method. 
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II 


PART I ALLY- PARABOLIC  PLOWS 


As  a  representative  example  of  high  Reynolds-number  flow  on 
a  three-dimensional  body,  we  consider  the  flow  around  a  ship 
hull.  Experiments  and  calculations  indicate  that  first-order 
boundary-layer  equations  adequately  describe  the  flow  over  a 
large  part  of  the  hull  if  phenomena  associated  with  the  free 
surface  are  excluded.  The  assumptions  embodied  in  boundary-layer 
theory  begin  to  breakdown  gradually  over  the  stern,  in  a  region 
which  is  typically  10  to  20  percent  of  the  ship  length.  The 
experimental  information  pertaining  to  the  evolution  of  the  flow 
over  the  stern  and  in  the  near  wake  has  been  reviewed  by  Patel 
(1982).  Among  the  features  which  differentiate  the  stern  flow 
from  the  boundary-layer  flow  over  the  hull  are  (a)  a  rapid  thick¬ 
ening  of  the  viscous  layer,  (b)  variation  of  pressure  across  the 
layer,  implying  a  strong  viscous- i nv i sc  id  interaction,  (c)  the 
development  of  a  large  longitudinal  vorticity  component  which  may 
or  may  not  lead  to  a  free-vortex  type  of  separation,  and  (d)  a 
general  reduction  in  the  level  of  turbulence.  The  wake  of  the 
ship  evolves  gradually  from  this  flow  and  also  shows  basically 
the  same  features.  An  important  observation  that  also  emerges 
from  the  experiments  is  that  there  is  usually  no  region  of  flow 
reversal  in  the  direction  of  ship  motion.  The  flow  features 
noted  above  are  precisely  those  which  can  be  addressed  by  the 
partially-parabol ic  equations.  In  particular,  these  equations 
can  be  used  to  describe  the  flow  in  the  region  between  the  thin 
boundary  layer  upstream  and  the  wake  far  downstream  of  the  hull. 

Extensive  regions  of  thick  boundary  layers,  strong  viscous- 
inviscid  interaction  and  longitudinal  vortex  formation  are  also 
observed  on  bodies  of  revolution  at  incidence  (see,  for  example, 
Patel  and  Baek,  1983).  For  boat-tailed  bodies  at  small  to  moder¬ 
ate  incidences,  longitudinal  flow  reversal  may  be  absent  and 
again  the  part ial ly-parabol ic  approximations  can  be  used  to  des¬ 
cribe  the  development  of  the  vortices  on  the  body  and  in  the 
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In  axisymmetric  flow,  the  longitudinal  vortices  are  of 
course  absent,  but  measurements  in  the  thick  boundary  layer  over 
the  tail  and  in  the  near  wake  of  five  different  bodies  (Patel  et 
al.,  1974;  Patel  and  Lee,  1977;  Huang  et  al.,  1979,  1980)  indi¬ 
cate  all  of  the  remaining  features  of  the  flow  over  a  ship 
stern.  As  in  the  three-dimensional  examples,  the  boundary  layer 
thickens  as  a  result  of  the  convergence  of  the  streamlines  in 
planes  parallel  to  the  surface  (or  due  to  the  decrease  in  the 
cross-sectional  area  of  the  body)  and  leads  to  strong  viscous- 
inviscid  interaction.  Here  again,  the  part ial ly-parabol  ic  (or 
the  identical  thin-layer)  equations  can  be  employed  to  extend  the 
boundary-layer  solutions  into  the  wake.  Many  different  methods 
have  been  developed  for  this  simple  class  of  flows  (see  Nakayama 
et  al.,  1976;  Dyne,  1978;  Geller,  1979;  Dietz,  1980;  Hoffman, 
1980,  1982;  Huang  et  al.,  1980;  Muraoka,  1980;  Zhou,  1982;  Hogan, 
1983;  Markatos,  1984).  All  except  the  last  four  of  these  employ 
the  classical  iterative  approach  to  account  for  the  viscous- 
inviscid  interaction. 

Two-dimensional  external  flows  in  which  the  part ial ly-para- 
boLic  approximations  mav  be  applied  are  those  involving  a  con¬ 
fluence  of  two  or  more  shear  layers.  Practical  examples  are  the 
flow  at  the  trailing  edge  and  in  the  wake  of  an  airfoil  in  the 
absence  of  separation,  the  interaction  between  wakes  and  boundary 
layers  as  on  multi-component  airfoils  and  downstream  of  wing-body 
junctions,  and  the  initial  regions  of  jets,  wall  jets  and  mixing 
layers.  The  near  wake  of  a  flat  plate  represents  the  limiting 
case  and  the  simplest  example  insofar  as  the  viscous- i nv i sc  id 
interaction  is  weak  and  confined  to  a  very  small  region  close  to 
the  trailing  edge.  Two-dimensional  tra i 1 i ng-edge  and  near-wake 
flows  have  been  treated  by  a  variety  of  methods,  including 
numerical  solutions  of  the  full  Nav ier-Stokes  or  Reynolds  equa¬ 
tions  (Hah  and  Lakshminarayana ,  1982;  Horstman,  1983;  Adair  et 
al.,  1983),  solutions  of  the  pa r t i a  1 ly-pa r abo 1 i c  equations  (Rubin 
and  Reddy,  1983),  and  methods  which  couple  the  external  inviscid 


flow  to  the  boundary  layer  and  wake  flow  through  successive  iter¬ 
ations  ( Le  Balleur,  1983;  Bradshaw  et  al.,  1983;  Vatsa  and 
Verdon,  1983).  McDonald  and  Briley  (1983)  have  reviewed  the 
recent  developments  in  the  latter  area  and  the  extensions  of  such 
techniques  to  treat  flows  with  small  regions  of  separation  using 
the  so-called  inverse  boundary-layer  calculations. 

The  method  to  be  presented  here  is  concerned  with  the  solu¬ 
tion  of  the  partially-parabolic  Reynolds-averaged  equations  for 
three-dimensional  flows  of  an  incompressible  fluid.  Although  our 
primary  interest  lies  in  the  flow  over  ship  sterns,  with  all  the 
complicating  factors  noted  above,  the  method  will  be  used  also 
for  the  solution  of  some  two-dimensional  and  axisymmetric  flows 
in  order  to  evaluate  its  detailed  performance. 

Previous  solutions  of  the  partially-parabolic  equations  in 
three-dimensional  external  flows  are  quite  limited.  Abdelmeguid 
et  al.  (1979)  presented  the  first  application  to  ship  hulls  and 
Markatos  et  al.  (1980),  Muraoka  (1980,  1982),  Tzabiras  and  Louka- 
kis  (  1983)  and  Tzabiras  (  1983)  have  pursued  it  further.  As  is 
the  case  with  many  other  applications  in  internal  flows,  all  of 
these  use  essentially  the  same  numerical  scheme,  based  on  the 
work  of  Patankar  and  Spalding  (1972),  and  the  k-e  model  of  turbu¬ 
lence  with  standard  wall  functions.  A  brief  review  of  these 
references  is  useful  to  indicate  the  present  status  of  the  sub¬ 
ject  and  also  to  point  out  the  major  features  of  the  new  method. 


II. 1  The  Partially-Parabolic  Approximations 

Strictly  speaking,  the  partially-parabol  ic  approximations 
are  physically  most  realistic  only  if  the  diffusion  in  the  direc¬ 
tion  of  the  local  mean  streamline  is  neglected.  In  all  previous 
calculations,  however,  the  part ial ly-parabol ic  equations  have 
been  derived  by  neglecting  the  stress  gradients  in  a  primary, 
oredomi nant-f low  direction  which  is  chosen  a  priori,  usually  the 
free-stream  direction  or  the  direction  of  ship  motion. 
Obviously,  significant  errors  may  arise  in  regions  where  the 


local  flow  is  not  in  the  primary  direction,  e.g.  close  to  a 
curved  surface  or  where  the  secondary  motion  is  not  weak.  If  the 
choice  of  the  primary  direction  is  postponed  until  a  suitable 
coordinate  system  is  selected,  some  of  the  terms  implicitly  neg¬ 
lected  in  the  previous  work  can  be  retained  without  altering  the 
part ial ly-parabo  1  ic  nature  of  the  equations.  This  is  discussed 
in  Section  1 1  I. 3. 

11.2  The  Coordinate  System 

In  all  previous  solutions  for  ship  hul Is,  the  Reynolds  equa¬ 
tions  are  first  written  in  the  cylindrical  polar  coordinates 
(x,r,0),  with  x  in  the  direction  of  ship  motion.  After  neglect¬ 
ing  the  x-derivatives  of  thr  stresses,  these  are  then  transformed 
such  that  the  hull  surface  becomes  one  of  the  coordinate  sur¬ 
faces.  Abdelmeguid  and  Muraoxo  employed  analytical  transforma¬ 
tions  which  led  to  a  nonorthog  via  1  distorted  polar  coordinate 
system  in  the  (r,0)  plane,  while  Tzabiras  and  Loukakis  employed 
conformal-transformat  in.  techniques  to  generate  orthogonal  co¬ 
ordinates  in  this  plane.  , n  both  cases,  the  coordinate  lines  in 
the  (x,r)  planes  are  r  orthogonal.  The  improved  results  ob¬ 
tained  by  the  latt-.  authors  have  been  attributed  to  this 
difference  in  the  coord  i  na  t.es  . 

Recently,  r'hen.;  anil  Patel  (L983)  have  suggested  the  use  of 
name  r  i  ca  11  '/-goner  a  -  ed  body-fitted  coordinates  for  ship  flow  cal- 
eul at  ions.  Their  method  has  been  extended  in  the  present  work  to 
a  more  general  class  of  coordinate  systems,  and  the  numerical 
procedure  for  the  calculation  of  the  coordinates  has  been  sim¬ 
plified  and  made  more  efficient.  The  basic  theory  is  discussed 
in  Section  IT  1.1  and  some  examples  are  presented  in  Section 
IV.  I  .  This  general  and  flexible  grid-generation  method  facili¬ 
tates  the  evaluation  of  the  influence  of  coordinate  systems  on 
the  numerical  solutions  and  also  enables  the  method  to  be  em¬ 
ployed  for  other  shapes. 


II. 3  Turbulence  Model 


As  in  the  previous  studies,  the  present  method  employs  the 
standard  two-equation  k-e  turbulence  model  and  uses  wall  func¬ 
tions  to  avoid  integration  of  the  equations  through  the  sub¬ 
layer.  One  of  the  objectives  of  the  research  is  to  evaluate  the 
efficacy  of  this  method,  and  the  role  of  the  turbulence  model,  in 
general,  in  the  flow  over  trailing  edges  and  ship  sterns.  This 
can  be  accomplished  by  detailed  comparisons  with  experimental 
data  only  after  confirming  the  accuracy  and  reliability  of  the 
numerical  solution  procedures.  Although  the  turbulence  model 
equations  are  given  in  Section  III,  specific  comments  concerning 
the  boundary  conditions  and  the  evaluation  of  the  model  are  post¬ 
poned  to  the  discussion  of  the  results  for  particular  cases. 

II. 4  Solution  Domain  and  Boundary  Conditions 

In  practical  applications  of  the  part ial ly-parabo 1 ic  equa¬ 
tions,  the  upstream  boundary  of  the  solution  domain  is  placed  at 
a  station  where  the  boundary  layer  is  thin  and  can  be  calculated 
by  first-order  theory.  The  location  of  downstream  boundary 
should  be  "far  enough"  from  the  body  such  that  the  local  flow  is 
parabolic,  i.e.  no  influence  propagates  upstream  from  there. 
Some  of  the  earlier  solutions  were  obtained  with  the  downstream 
boundaries  located  in  the  near  wake,  and  therefore  contain  an 
unknown  influence.  The  coordinates  and  the  numerical  methods 
used  here  enable  the  downstream  boundary  to  be  placed  several 
body-lengths  downstream  without  a  large  penalty  in  computer  time. 

The  extent  of  the  solution  domain  normal  to  the  body  is 
usually  taken  to  be  of  the  order  of  the  thickness  of  the  viscous 
flow  and  the  conditions  at  this  outer  boundary  are  specified  from 
poten t i al- f low  solutions.  In  view  of  the  interaction  between  the 
viscous  and  invisid  flows,  several  iterations  are  required  to 
obtain  a  converged  solution.  While  this  procedure  is  commonly 
employed  in  two-dimensional  and  axisymmetric  flows,  previous 
solutions  for  three-dimensional  flows  have  been  obtained  without 
any  iteration,  i.e.  by  simply  specifying  the  external  conditions 


from  the  potential  flow.  Although  usual  methods  of  displacement- 
thickness  or  equivalent  surface  singularities  can  be  used  to 
account  for  the  influence  of  the  viscous  flow  on  the  external 
flow,  the  computation  times  have  thus  far  proved  to  be  quite 
prohibitive  for  fully  three-dimensional  interactive  solutions. 

One  of  the  important  features  of  the  present  method  is  that 
the  external  boundary  of  the  computation  domain  is  placed  at  a 
large  enough  distance  from  the  body  so  that  uniform-flow  condi¬ 
tions  can  be  imposed  there.  This  captures  the  entire  zone  of 
v iscous- i nv  i  sc  id  interaction  and  avr  ids  the  need  for  potential- 
flow  and  iterative  solutions  altogether.  As  we  shall  see  from 
the  examples  presented  later  on,  this  particular  feature  is  made 
possible  by  coordinate  stretching  and  the  properties  of  the  "fi¬ 
nite-analytic"  numerical  method  that  is  employed. 

II. 5  Numerical  Scheme 

For  two-dimensional  and  axisymmetric  problems,  different 
numerical  schemes  have  been  used  to  solve  the  partially-parabol ic 
equations.  The  methods  also  differ  in  the  manner  in  which  the 
pressure  field  is  calculated  to  allow  for  the  upstream  propaga¬ 
tion  of  information.  This  has  usually  involved  several  sweeps 
through  the  computation  domain  marching  from  the  upstream  to  the 
downstream  boundaries.  In  the  few  applications  to  three-dimen¬ 
sional  external  flows  mentioned  earlier,  the  implicit  finite- 
volume  numerical  scheme  of  Patankar  and  Spalding  (1972)  has  been 
used  together  with  the  SIMPLE  algorithm  of  Patankar  (1980)  for 
pressure  correction  due  to  the  read  •  availability  of  commercial 
computer  programs.  Tn  these  calculations,  an  initial  pressure 
field  is  obtained  from  the  potential-flow  solution  and  the  pres¬ 
sure  is  corrected  in  successive  sweeps  of  the  viscous  solution. 
Pratap  (197b)  and  others  have  noted  a  rather  slow  rate  of  conver¬ 
gence.  Use  of  this  pressure-correction  procedure  in  our  preli¬ 
minary  studios  led  to  an  oscillatory  pressure  field  which  lasted 
for  several  hundred  sweeps.  rt  was  therefore  found  necessary  to 
adopt  alternative  methods. 


The  present  numerical  scheme  has  three  distinctive  fea¬ 
tures.  First,  it  employs  the  recently-developed  finite-analytic 
numerical  scheme  (Chen  et  al.,  1982)  which  is  quite  different 
from  the  conventional  finite-difference  and  finite-element  me¬ 
thods.  In  fact,  it  can  be  viewed  as  a  combination  of  the  two 
approaches.  It  uses  analytical  solutions  of  the  local ly- 1 i near- 
ized  convective  transport  equations  for  each  numerical  cell  in 
the  computation  domain  and  combines  them  to  obtain  the  global 
solution  by  conventional  methods.  Secondly,  the  pressure  correc¬ 
tion  is  obtained  by  a  two-step  procedure  which  uses  a  modified 
version  of  the  SIMPLER  algorithm.  This  improves  the  convergence 
properties  of  the  overall  solution.  Thirdly,  the  method  is  for¬ 
mulated  for  the  more  general  unsteady  three-dimensional 
problem.  While  this  leads  to  a  small  increase  in  computer- 
storage  requirements,  it  facilitates  the  management  of  informa¬ 
tion  and  enables  future  extensions  to  unsteady  and  f ul ly-el 1 ipt ic 
flows.  For  steady  flows,  the  converged  solution  is  also  obtained 
by  time  marching  but,  in  this  case,  time  serves  as  an  iteration 
parameter.  Details  of  the  numerical  scheme  are  discussed  in 
Section  IV. 
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EQUATIONS  AND  COORDINATE  SYSTEMS 


Since  most  of  the  applications  presented  in  this  report  are 
associated  with  external  flows  past  ax i symmetric  and  three-dimen¬ 
sional  bodies,  such  as  ship  forms,  it  is  convenient  to  choose 
cylindrical  polar  coordinates  as  the  basic  coordinate  system  in 
the  physical  domain.  Formulations  start  in-;  with  other  basic 
coordinate  systems  can  be  derived  in  a  similar  manner  (see  Appen¬ 
dix  I  )  . 

Consider  the  equations  of  mot i  n  in  cylindrical  coordinates 
(x,r,9)  for  unsteady,  three-  1  imens ;  >nal,  incompressible  flows. 
The  exact  Rey nolds-ave raqed  equations  °f  continuity  and  momentum 
of  the  mean  flow  in  dimension!  ,;s  form  are. 
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where  x  =  — ,  r  =  — ,  0  are  the  dimensionless  coordinates  normal- 

L  L 

ized  by  a  characteristic  length  L,  and  t  is  the  time  normalized 
by  L/Uq.  U,  V,  W  are,  respectively,  the  longitudinal,  radial  and 
circumferential  components  of  mean  velocity  normalized  by  the 
characteristic  ^vglocity  UQ.  p  is  the  pressure  normalized 
by  pU0^*  Re  =  ~ is  the  Reynolds  number  defined  in  terms  of  UQ , 
L  and  molecular  kinematic  viscosity  v.  The  barred  quantities 
uu  ,  uv ,  etc.  are  the  Reynolds  stresses,  normalized  by  u^. 

In  the  present  study,  the  two-equation  k-e  turbulence  model 
is  used  to  model  the  Reynolds  stresses.  Each  Reynolds  stress  is 
related  to  the  corresponding  mean  rate  of  strain  by  an  isotropic 
eddy  viscosity  v  as  follows: 
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The  eddy  viscosity  is  related  to  the  dimensionless  turbulent 
kinetic  energy  k,  and  its  rate  of  dissipation  e,  by 


v.  =  c  — 
t  ye 


(  I  I  1-6  ) 


where  c  (=0.09)  is  a  constant,  and  k  and  c  are  governed  by  the 
convective  transport  equations 


+  u 


+  V 


at 


a  x 


ar 


+  7  M 


l  a 
+  — 


1 


Tx 

i  a 


o  r  _ 
k  e  f  f 


75^ 


1 


ak  , 


ak  . 

r  ar  '  o ,  K  , ,  r  3  r  +  2  38  'O  R  36 

k  eff  r  k  e  f  f 


1  +  G  -  E 


3  t 

TF 


a  t 
a  x 


l  a 


a  l  w 

a  e 

a 

r  .  i  ^ 

-  +  — 

j  r  r 

7e 

ax 

v  o  R  3x  > 

c  eff 

i 

a  e  . 

l  a  f  l 

a  e  ^ 

0  R  _  £ 
c  elf 

IT  n 

d  r 

+ 

2  3  0  v  o  R  ,, 
r  e  ef  f 

3  0  ; 

(  i  i  g 

t  1  k 


C  t 2  k 


where  G  is  t. he  turbulence  generation  term: 
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The  effective  Reynolds  number,  R^^,  is  defined  as 
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The  constants  in  these  equations  were  taken  as  C  =  0.09,  C  = 

pel 

1.44,  C  =  1.92,  0  =  1.0,  and  o  =  1.3.  We  note  again  that 
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v  ,  k,  and  f-  are  all  made  dimensionless  by  the  characteristic 
velocity  and  length  scales  U0  and  L,  respectively. 

Using  equations  (III-5,  I  1 1  —  6  and  III  — 10) ,  the  momentum 
equations  (II  1-2)  through  ( 1 1  I  —  4 )  and  the  turbulence-model  equa¬ 
tions  (  I  I  I  —  7 )  and  (  1 1 1  —  8 )  are  transformed  to 
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It  is  convenient  to  rewrite  the  governing  transport  equations 
(  III— 11)  through  (  III  — 15)  in  the  following  general  form 
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where  $  may  represent  any  one  of  the  convective  transport  quanti¬ 
ties:  U,  V,  W,  k  or  t,  and  the  subscripts  x,r,0  denote  deriva¬ 
tives.  The  corresponding  coefficients  a  ,  b  ,  c  and  d  are 
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Rquations  ( 1 1 r  —  1 1 )  through  < III  — 15)  are  coupled,  nonlinear,  par¬ 
tial  differential  equations  and,  together  with  the  continuity 
equation  (III-l),  are  sufficient,  in  principle,  to  solve  for  the 
six  unknowns  p,  U,  V,  W,  k  and  e  when  proper  initial  and  boundary 
conditions  are  specified.  It  should  be  noted  that  these  are 


still  the  complete  elliptic  equations  and  no  approximations  have 
been  made  other  than  those  inherent  in  the  turbulence  model. 

ITI.l  Body-Fitted  Coordinates 

For  most  practical  applications,  the  geometry  of  the  body  is 
usually  quite  complex  and  the  cylindrical  or  other  basic  ortho¬ 
gonal  coordinates  are  not  the  most  convenient  ones  to  use.  It  is 
therefore  desirable  to  introduce  analytic  or  numerical  coordinate 
transformations,  which  facilitate  the  applications  of  the  boun¬ 
dary  conditions  and  also  simplify  the  computational  domain  in  the 
transformed  plane. 

As  mentioned  earlier,  two  types  of  coordinate  systems,  name¬ 
ly  a  distorted  polar  coordinate  system  obtained  by  analytic 
transformation  and  that  based  on  conformal-transformation  at  each 
section,  have  been  employed  in  previous  investigations  associated 
with  ship  hulls.  From  the  limited  calculations  performed  with 
these,  it  appears  that  the  solutions  are  greatly  dependent  upon 
the  choice  of  coordinates  since  all  have  used  essentially  the 
same  equations,  turbulence  model  and  numerical  method.  There¬ 
fore,  it  is  important  to  investigate  other  coordinates,  and  tech¬ 
niques  for  coordinate  generation,  so  that  accurate  and  efficient 
flow  calculations  can  be  performed.  Here  we  shall  follow  the 
increasingly  popular  method  of  numerical  grid  generation,  as 
suggested  recently  by  Cheng  and  Patel  (IS) 83),  since  it  offers  the 
advantages  of  generality  and  flexibility  and,  most  importantly, 
transforms  the  computation  domain  into  a  simple  rectangular  re¬ 
gion  with  equal  grid  spacing. 

In  the  numerical  grid-generation  technique,  we  seek  a  co¬ 
ordinate  system  for  the  numerical  analysis  of  the  flow  in  the 
domain  D  shown  in  Fig.  1.  This  domain  is  bounded  by  an  arbitrary 
hull  surface  S,  the  ship  centerplane  C,  the  free  surface  or  water 
plane  W,  the  upstream  and  downstream  sections  A  and  R,  respec¬ 
tively,  and  an  external  boundary  I.  Section  A  may  be  located  at 
a  hull  section  where  the  boundary  layer  is  thin,  and  B  may  he 
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placed  at  a  section  in  the  far  wake.  The  choice  of  the  external 
boundary  E  is  arbitrary.  The  basic  idea  of  a  boundary-conforming 
curvilinear  system  is  to  find  a  transformation  such  that  the 
curvilinear  boundary  surfaces  of  the  physical  domain  D  in  cylin¬ 
drical  or  in  any  other  basic  orthogonal  coordinate  system,  say 
1  2  i 

(x  ,  x  ,  x  ),  are  transformed  into  boundaries  of  a  simple  rec¬ 
tangular  domain  in  the  computational  space  (£,  n ,  c)  shown  in 
Fig.  2. 

With  the  values  of  the  curvilinear  coordinates  specified  on 
the  boundaries  of  D,  it  then  remains  to  generate  the  values  of 
these  coordinates  in  the  interior  of  D.  This  is  a  boundary-value 
problem  in  the  physical  field  with  the  curvilinear  coordinates 

(4,  n ,  O  as  dependent  variables  and  the  orthogonal  coordinates 
12  3 

(x  ,x  ,x  )  as  the  independent  variables,  with  boundary  conditions 
specified  on  the  curved  boundaries.  Thus,  a  system  of  elliptic 
partial  differential  equations  can  be  used  to  generate  the  co¬ 
ordinates  since  the  field  solution  of  such  a  system  is  determined 
entirely  by  the  boundary  conditions.  However,  the  elliptic  sys¬ 
tem  must  be  chosen  such  that  it  precludes  the  occurrence  of 
extrema  in  the  interior  of  the  domain  and  assures  a  one-to-one 
napping  between  the  physical  and  the  transformed  planes. 

For  the  general  three-dimensional  but  simply-connected  do¬ 
main  of  interest  here,  a  set  of  Poisson  equations  of  the  form 

=  fl(S,  n,  c) 

v2n  =  f2U,  n  ,  a  (  1 1 1-19  ) 

v2c  =  f3(C,  n ,  O 

or  simply 

V2  C 1  =  f1,  i  =  1,  2,  3  (HI-20) 

with  =  f,,  f,2  =  n,  f,3  =  c 
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may  be  taken  as  the  basic  coordinate  generating  system.  Here, 
V  is  the  Laplacian  operator  in  orthogonal  coordinates  x1.  The 
nonhomogeneous  source  functions  f1  may  be  assigned  appropriate 
values  to  yield  the  desired  concentration  of  coordinate  sur¬ 
faces.  The  choice  of  these  functions  for  specific  applications 
will  be  discussed  later  in  Section  V.  Equation  (III-20)  is 

subject  to  either  Dirichlet  or  Neumann  boundary  conditions  on  the 
boundary  surfaces,  which  are  surfaces  of  constant  £ 1 . 

Since  it  is  very  desirable  to  perform  all  numerical  computa¬ 
tions  in  the  transformed  (£,  n,  5)  plane  with  equal  grid  spacing, 
i.e.  A  £  =  An  =  At;  =  1,  equation  (III-20)  is  cumbersome  to  use. 
It  is  more  convenient  to  invert  it  and  solve  for  the  orthogonal 
coordinates.  In  other  words,  the  dependent  and  independent 

variables  are  interchanged  so  that  the  orthogonal  coordinates 
12  3 

(x  ,x',x  )  in  the  physical  plane  become  the  dependent  variables, 
with  the  curvilinear  coordinates  (£,  n,  c)  as  the  independent 
variables.  The  bounda ry- va 1 ue  problem  in  the  transformed  field 
then  involves  generating  the  values  of  the  orthogonal  coordinates 
x 1  -  xl  (£,  n,  ?)  in  the  interior  from  the  specified  boundary 

values  of  xl  on  the  rectangular  boundary  surfaces  of  the  trans¬ 
formed  field.  Since  the  boundaries  in  the  transformed  plane  are 
all  rectangular  (constant  £,  n,  or  £  plane),  these  computations 
are  carried  out  on  a  cubic  grid  regardless  of  the  shape  of  the 
physical  boundaries. 

In  obtaining  the  inverse  transformation  of  equation  (III- 

12  3 

20),  several  general  relations  between  the  physical  (x  ,  x',  x  ) 
and  the  transformed  (  £  ,  n,  r.)  coordinates  are  required.  The 
basic  expressions  may  be  found  in  many  reference  books,  for 
example  Aris  (1971)  and  Lass  (1975),  although  some  of  the  rela¬ 
tions  are  not  explicitly  given.  For  the  sake  of  completeness, 
therefore,  some  important  relations  are  given  in  Appendix  I  with 
specific  reference  to  the  transformation  between  any  orthogonal 
coordinates  ( x  ^  ,  x^,  x^)  and  general  coordinates  (  £  ^  ,  £^f  £^)  - 
(£,  n,  £).  The  latter  are  not  necessarily  orthogonal.  With 
these  transformation  formulae,  equations  ( III  —  20)  become: 
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h^  are  the  metric  coefficients  in  the  chosen  orthogonal  coordi¬ 
nates  x1,  and  g1-'  is  the  conjugate  metric  tensor  in  the  trans¬ 
formed  coordinates  4  1  (=  n,  C )  .  Note  that,  for  the  cylindrical 
polar  coordinates  (x,  r,  9)  chosen  here,  the  metrices  h|  are  (1, 
1,  r),  respectively. 

Equations  (III  —  21)  can  be  solved  numerically  in  the  trans¬ 
formed  domain  (S,  n,  C)  when  proper  boundary  conditions  are  spec¬ 
ified  on  all  boundary  surfaces  (i.e.  constant  £,  n  and  c  ).  If 
12  3 

f  =  f  =  f  =0,  the  transformation  is  said  to  be 
homeomorph ic .  In  general,  however,  non-zero  values  may  be  as¬ 
signed  to  these  functions  to  exercise  control  over  the  grid  dis¬ 
tribution.  Solutions  of  Equations  ( III  —  21 )  to  obtain  numerical 
grids  for  particular  shapes  are  presented  in  Section  V. 


III. 2  The  Reynolds  Equations  in  Transformed  Coordinates 

In  the  previous  section,  the  numerical  grid-generation  tech¬ 
nique  based  on  a  coordinate  transformation  from  an  orthogonal 
curvilinear  system  to  a  general  curvilinear  system  was  out¬ 
lined.  It  then  remains  to  transform  the  governing  equations  of 

motion  to  the  transformed  plane.  Although  it  is  possible  to 

1  O  1 

transform  both  the  independent  (x  ,  x  ,  x  ,  t)  and  the  dependent 
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(U,  V,  W,  k,  e)  variables  to  the  ( £  ,  n,  ?)  coordinates,  we  will 
consider  only  the  tranforraat ion  of  the  independent  variables, 
leaving  the  velocity  components  U,  V,  W  in  the  original  (x  ,  x  , 
x^)  coordinates  in  the  physical  plane.  This  approach  makes  it 
possible  to  retain  a  certain  amount  of  physical  intuition  in  the 
approximations  to  be  made  later. 

As  mentioned  earlier,  the  cylindrical  polar  coordinates  (x, 
r,  9)  appear  to  be  most  convenient  for  the  description  of  the 
flow  field  around  axisymmetric  bodies  and  ship  forms.  Thus,  in 
this  study,  equations  (  1 1 1  —  1 )  and  ( III  — 11)  through  ( III  — 15)  will 
be  used  as  the  basic  equations  in  the  physical  plane  to  derive 
the  equations  in  the  transformed  domain  (£,  n,  C).  Transformed 
equations  resulting  from  other  basic  orthogonal  coordinates, 
which  may  be  useful  in  other  applications,  are  given  in  Appendix 
I. 

Using  the  general  expression  of  the  gradient,  divergence  and 
Laplacian  given  in  Appendix  I,  equations  (TIT  — 11)  through  (  [  r  I  — 
15)  for  unsteady  three-dimensional  flow  can  be  transformed  to  the 
following  form 
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and  a  .  h  c  and  dA  are  as  defined  in  equation  (II 1-17).  Here- 
4>  4>  't'  t 

after,  the  subscripts  (  f, ,  n ,  C )  on  t  (=[J,V,W,k,  e)  and  denote 

derivatives. 

liquation  (ITl-22)  is  identical  with  equation  (A-84)  in 
Appendix  I.  It  can  be  rearranged  into  a  general  convect ive- tran¬ 
sport  equation  of  the  form 
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and  the  source  functions  are  defined  in  equations  (A-85a) 


through  (A-85f)  in  the  transformed  plane.  For  the  cylindrical 


polar  coordinates  (x,  r,  9)  chosen  here,  the  only  nonzero  k 


and  a ..are  k 


i  3 
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functions  are  greatly  simplified: 
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(  III-26 a) 


eff  {-  f~  +  J  (b2P^+  b2Pn  +  b2Pc>  +  fj  (b2kC+  b2kn+  b2V 
J2  (blvt,C+  bl  Ut,n+  blvt,s } (b2US+  b2Un+  b2U?} 

-  J  (b3  vt,  C+  b3  vt ,  n+  b3  vt,C^J  (b2V  b2Wn 
+  b2  VV  ~  ?]>  +  TJ  <b3w£+  b3Wri+  b3w?)  +  ^2  (III-26h) 


eff{  r  +  J  (b3p4+  b3pn+  b3  PC ]  +  3J  (b3k£+  b3  kn+  b3k;} 

'  T2  <bK,C+  bl  vt ,  n+  b?\,C)(b3V  b3b  n+  b3'V 

yj 

-  -7  ( bi  v  ,+  b2v.  +  b^v.  W-7  ( b^Vr+  b^V  +  bi?v  )  -  -] 

-  M  <b3Vt,t+  bK,n  +  b3“t,C)l  -  73  <b3VG  b3Vn+  b3VC> 


(  1 1 1-2fic ) 


o  R  (G  -  e: 
k  eff 


(  1 1 1-26<1 ) 


a  R  ^  (C_,  t  G  -  C. 


(  1 1 1-26e ) 


G  =  vt  {72  (blV  W  biU,)2+  72  (b2Vr/  b2Vn+  b2V  ' 

J  iJ 

+  2  [j  (b*W^  b2W,+  b*W;>  +  |]2 
+  7  <blV  blVn  +  blV  b2U5  +  b2U„  +  b2Ut»2 

J 

+  ~  ^blWS+  biwn+  blwc+  b3U£  +  b3Un  +  b3Uc' 
lJ 

^  jL  2  ^  |  2  ^  ri  2 

+  [j  (b2W^+  b2Wn  +  t>2W ,  +  b3V^+  b3Vn+  b3V;)  -  -]  }  (  1 1 1- 26  f  ) 

The  radius  r,  the  geometric  coefficents  b^  and  g13,  and  the 

Jacobian  J  which  appear  in  the  above  equations  are  functions  of 

the  coordinates  only.  When  either  analytic  or  numerical 

transformations  are  employed  to  generate  the  grid  distributions, 

r  =  r(£,  n,  5)  is  known  in  the  transformed  plane.  The  geometric 

coefficients  b^  and  g13  for  the  present  transformation  are 

1  9 

obtained  from  equations  (A-35),  (A-40),  and  (A-42)  with  (x  ,  x  , 
x3)  =  (x,  r,  0)  and  ( h 1 ,  h2,  h3)  =  (1,  1,  r),  i.e. 
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is  the  determinant  of  the  metric  tensor  g^j.  Also,  the  Jacobian 
can  be  expressed  as 


r  9 ,  r0„  r  9  r 
C  n  S 


( III-31  ) 


It  should  be  remarked  here  that  equation  (III. 24)  is  still 
fully  elliptic  in  space  and  (U,  V,  W)  represent  the  velocity 
components  in  the  physical  cylindrical  polar  coordinates 
(x,  r,  9).  Equations  ( III  —  24)  through  (III-26),  together  with 
the  equation  of  continuity,  i.e. 
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(  III-32) 


(see  equation  A-66)  are  the  Reynolds-averaged  Nav ier-Stokes 
equations  for  unsteady,  three-dimensional  flows. 

Before  considering  the  part ia 1 ly-pa rabol i c  approximations  it 
is  useful  to  write  down  the  full  Reynolds  equations  for  two 
special  cases,  namely  axisymmetric  flow,  which  illustrates  the 
simplifications  that  arise,  and  three-dimensional  flow  in  the 
distorted  cylindrical  polar  coordinate  system  which  has  been  used 
by  some  previous  investigators. 


(  i )  Axisymmetric  flow 

In  this  case,  U  =  V  =  k  =  e  =p=x=r=  9  =  0  =  0  and 

cccc^cccS  n 

0^  =  1.  Therefore,  equation  (III-24)  reduces  to 
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The  coefficients  b^  and  glJ  in  the  above  equations  can  be 
obtained  easily  from  equations  (III-27)  through  (III  —  31 )  by 


noting  that  9^  =  9  0  and  6^  =  1,  and  are  therefore  given  by 
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( III- 36) 
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and  the  Jacobian  is 


J  =  •'g  *  r(xrr„"  x  rr) 


(  I  1 1  -  38 ) 


These  equations  are  used  later  to  calculate  the  flow  in  the 
thick  axisymmetric  boundary  layer  over  the  tail,  and  in  the  near 
wake,  of  bodies  of  revolution  for  which  extensive  data  are 
available. 

(ii)  Equations  in  Distorted  Cylindrical  Polar  Coordinates 


The  distorted  coordinates  employed  by  Abdelmeguid  et  al. 
(1979)  and  Muraoka  (1980,  1982)  correspond  to  the  analytical 

transformations 
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(  III-39) 


where  x  is  measured  along  the  direction  of  ship  motion,  rn  and  rs 
represent  the  radial  distances  to  the  preselected  outer  and  inner 
(hull)  boundaries,  respectively,  and  ^  =  n/2  is  the  extent  of  the 
domain  in  the  circumferential  direction,  i.e.,  from  the  waterline 
to  the  keel.  Thus,  n  and  ?  vary  between  0  and  1,  and  boundary 
conditions  are  applied  along  n  =  0,  1  and  C  =  0,  1. 

When  the  abbrev iat ions 
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l  o 

are  introduced,  the  coefficients  a^  and  b^  in  the  transformed 
plane  can  be  written,  by  equations  ( A—  3  3)  and  ( III  —  27) , 


(  II 1-41 ) 
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where  b0  is  the  cofactor  of  a.  ,  such  that  a^  = 


J6.  ,  and  the 
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Jacobian  J  is  defined  by 


J  =  det  (a^  )  =  h3A 


The  metric  tensor  g—  is  given  by 
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with 


g  =  det  (g^ j )  =  J2  =  h2  A2 


and  the  conjugate  metric  tensor  g1-'  is 
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Mso,  the  functions  f 1 ,  f2  and  f3  are  respectively. 
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Therefore,  the  governing  equations  in  these  coordinates  can 
be  written 
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When  the  part ial ly-parabol ic  assumptions,  discussed  later, 
are  made  by  neglecting  the  stress  gradients  along  the  direction 
of  ship  motion  x,  the  above  equations  reduce  to  those  employed  in 
Abdelmeguid  et  al.  (1979)  and  Muraoka  (1980,  1982). 


Thus,  the  general  equations  (III-24)  and  (III-32)  can  be 
used  with  analytic  as  well  as  numerical  coordinate 
transformations.  When  other  analytic  or  numerical  coordinate 
transformations  are  considered,  one  needs  only  to  replace  the 
three  sets  of  coefficients  b^,  g1^  and  f1.  The  subsequent 
solution  procedure  remains  exactly  the  same.  This  demonstrates 
the  generality  and  flexibility  of  the  present  approach  in  setting 
up  the  governing  equations  and  coordinates. 


III. 3  The  Partially-Parabolic  Equations 

As  discussed  in  Section  IT,  the  part ial ly-parabol ic 
approximations  are  based  on  the  physical  assumption  that 
longitudinal  diffusion  due  to  viscosity  and  turbulence  may  be 
neglected  in  a  certain  class  of  flows  at  high  Reynolds  numbers. 
Noting  that  in  most  flows  of  practical  interest,  and  in 
particular  on  ship  sterns,  the  direction  of  the  mean  streamlines 
may  vary  significantly  from  point  to  point,  the  most  appropriate 
coordinate  system  which  captures  the  spirit  of  the  partially- 
parabolic  approximations  is  the  intrinsic  or  streamline 
coordinate  system.  While  such  a  system  has  been  used  recently  in 
the  simpler  case  of  axisymmetric  flows  (Zhou,  1982;  Hogan,  1983), 
it  is  not  a  very  practical  one  for  three-dimensional  flows.  In 


most  previous  applications  of  part ia 1 ly-parabol i c  equations, 
diffusion  has  been  neglected  in  some  preselected  primary  or  pre¬ 
dominant  flow  direction,  e.g.  in  the  direction  of  ship  motion. 
Thus,  all  x-derivatives  of  the  viscous  and  turbulent  stresses  are 
dropped  in  the  Reynolds  equations  (III-2)  through  ( 1 1 1  —  4 )  prior 
to  a  coordinate  transformation.  Obviously,  significant  errors 
may  arise  in  regions  where  the  direction  of  the  local  flow 
deviates  appreciably  from  the  selected  primary-flow  direction. 

In  addition,  the  omission  of  $  is  not  sufficient  to  render  the 

xx 

equations  partial ly-parabol ic  in  the  transformed  domain  unless  a 
particular  type  of  coordinate  system,  in  which  S=£(x),  is 
selected.  This  can  be  seen  clearly  from  the  following.  Since 
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the  viscous  diffusion  terms  u?  <t>  become 
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'=  xc*  y?+  v  etc* 


( ni-54) 


Therefore,  the  omission  of  does  not  strictly  imply  that  g^  = 

0.  In  order  to  render  the  transformed  equations  partially- 

paraholic,  one  must  select  S=S(x),  i.e.,  £  =  £  =  0,  so  that 

y  z 

there  is  no  contribution  of  4>  and  <j>  to  4>  _  _ .  In  other  words, 

yy  zz  ££ 

additional  approximations  must  be  made,  explicitly  or  implicitly, 
to  render  the  transformed  equations  par t ia 1 1 y-parabol ic  if  a 
general  coordinate  transformation  £=£(x,y,z)  is  to  be  used. 

In  the  present  work,  an  alternative  approach  is  employed  by 
introducing  the  part i a  1 ly-parabol ic  approximations  after 
transforming  the  f ully-e i 1 iptic  Reynolds  equations  into  the  body- 
fitted  coordinates.  Mathematically,  equations  ( 1 1 1  —  2  4 )  can  be 
rendered  pa r t i a  1 1 y-pa rabo 1 ic  by  simply  neglecting  the  terms 
involving  the  second  derivatives  with  respect  to  £,  the  marching 
direction,  to  obtain 


22,  3  3 

rJ  *nn  +  y  * 


tr  2  VP  2VP  VP  VP  s, 


( II 1-55) 


Where,  all  coefficients  remain  the  same  as  defined  earlier.  Note 
that  all  first  derivatives  in  the  £-  direction  arising  out  of  the 
diffusive  terms  and  coordinate  transformations  are  retained  in 
these  coefficients.  These  approximations  do  not  strictly  imply 
that,  the  physical  diffusion  in  the  £-  direction  is  neglected,  as 
the  following  shows. 


From  Equation  (III-54),  we  note  that,  even  for  laminar 

flows,  ug^4>^  does  not  represent  the  physical  diffusion  in  the 

4-  direction  unless  g^=g^  =  0  and  f1  =  0.  The  former  requires  an 
orthogonal  coordinate  system  (£,n,;)  in  the  transformed  plane, 
while  the  latter  places  a  restriction  on  the  grid  distribution. 
The  interpretation  of  <t>^  becomes  even  more  difficult  when 
turbulent  flows  are  considered.  For  example,  the  molecular  and 

turbulent  diffusion  of  x-momentum  in  the  x-direction  in  (x,y,z) 

or  (x,r,  0)  coordinates,  using  equation  (  III  —  5) ,  is  simply 

3  311  - ,  . .  .  32U  .  „  3ut  3U  2  3k 

Tt  t2v  7Z  ~  uul  -  2(  v  +  vt)  — +  2  jj-  37  ~  3  JZ 

3x 

In  our  general  notation,  this  is 


>  +  2v  4> 

XX  t ,X  X 


-  -r  k  ; 
3  x 


( 1 1 1- 56  ) 


Mathematically,  it  is  necessary  to  neglect  only  the  first  term  to 

render  the  x-momentum  equation  parabolic  in  the  x-direction.  The 

second  term  can  be  retained  and  treated  in  the  same  manner  as  the 

convective  term  u<t>  in  the  equations  while  the  third  can  be 

x 

absorbed  into  the  source  terms. 

Whether  the  present  approach  represents  an  improvement  over 
the  convectional  one  remains  to  be  determined  by  actual 
solutions.  Nevertheless,  two  important  features  of  this  approach 
are  noteworthy.  First,  invoking  the  partially-parabol ic 

approximations  in  the  body-fitted  coordinates,  in  which  £  is 
taken  along  the  body,  and  along  an  external  streamline  far  from 
the  body,  ensures  that  diffusion  is  neglected  in  a  direction 
roughly  coincident  with  the  local  streamlines  in  the  solution 
domain.  Secondly,  the  omission  of  the  fewest  terms  in  the 
complete  Reynolds  equations  facilitates  the  future  extension  of 
the  solution  procedure  to  the  f ul ly-e 1 1 i pt ic  equations.  In  fact, 
one  may  simply  retain  r  r  in  the  source  functions  S,  and  estimate 
their  contributions  from  the  solution  at  a  previous  time  step  or 


iteration.  This  is  then  a  f ul ly-el 1 ipt  ic  formulation  with  upwind 
differencing  for  first  derivatives  with  respect  to  £ . 


In  order  to  illustrate  the  physical  significance  of  the 
part ial ly-parabol ic  approximations  in  general,  it  is  instructive 
to  consider  the  following  linear  one-d imens ional  elliptic 
transport  equation  with  constant  positive  convective  coefficient 
A  and  source  function  g 


4> 


xx 


2  Ai>  +  g 
x  ^ 


subjected  to  the  boundary  conditions 


( 1 1 1-5  7 ) 


M-h)  =  i  ;  ❖  (  h )  = 


where  the  subscripts  U  and  D  denote  upstream  and  downstream, 
respect i vely .  The  analytic  solution  of  equation  ( III  —  57) ,  when 
evaluated  at  the  center  node  P  of  the  interval  -  h<x<h,  gives 
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(  1 1 1-58 ) 


This  indicates  that  the  influence  from  downstream  decreases 
rapidly  as  the  convective  coefficient  2A  (or  the  cell  Reynolds 
number)  becomes  large.  Furthermore,  the  derivative  at  point  P 
can  also  be  obtained  analytically  as 
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P 


2A(*  -  *  )  +  2gh 

2Ah  -  2 Ah 

e  -  e 


(  1 1 1-59 ) 


For  large  convective  coefficients,  this  simply  reduces  to 


♦x|p  -  -  2A  (111-60) 

in  agreement  with  equation  (It  1-57)  with  1  =  0. 

xx 

Thus,  with  the  present  pa rt i a  1 ly-pa rabo 1 i c  approximations 
(  $  ,  =  0),  the  s-  derivatives  associated  with  the  convective 


terms  ( <t>  ,  v  )  can  be  evaluated  using  backward  or  upwind 

S  t  f  S 

differences,  i.e., 


( 1 1 1-6 1 ) 


This  is  also  supported  by  the  analytic  solution  of  the  model 
equation  since,  from  equation  (III-58), 
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Ah  -Ah 
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2_ 

2A 


tanh  Ah 
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(  II 1-62  ) 


for  large  A,  in  agreement  with  the  asymptotic  behavior  noted 
above. 

Note  that  the  convective  coefficient  A  in  equations  (IIT-57) 
and  (  III  —  60)  may  contain  v  and  f^,  in  addition  to  U,  for 

t  fX 

turbulent  flows  in  general  coordinates.  Therefore,  the  second 
term  in  equation  (III-56)  should  be  retained  for  consistency 
although  4>  is  no  longer  the  physical  diffusion. 

With  the  £-  derivatives  in  the  general  transport  equations 
( III  —  55)  evaluated  by  upwind  differences,  it  is  possible  to  solve 
them  by  a  numerical  technique  marching  in  the  C-  direction, 
provided  the  pressure  gradients  in  the  source  terms  S^(4>  =  U,  V, 
W)  are  known.  The  pressure  is  obtained,  in  turn,  by  requiring 
the  velocity  field  to  satisfy  the  continuity  equation  ( III  —  32) . 
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IV.  NUMERICAL  SOLUTION  PROCEDURES 


In  this  section  we  consider  the  numerical  techniques  Cor 
grid  generation  and  the  solution  of  the  part ial ly-parabol  ic 
equations  (III-55) 

g2\n  +  =  2V?  +  2Vn  +  Vc  +  Vt  +  s<d  (IV~1) 

and  the  equation  of  continuity  (III-32) 

(bju  +  b^V  +  b^W)?  +  (b2U  +  b2V  +  b3W)n 

+  (b^U  +  b^V  +  b^W)  =  0  ( IV- 2 ) 

Details  of  the  numerical  method  used  in  grid  generation,  together 
with  some  examples,  are  presented  in  Section  IV. 1.  The  finite- 
analytic  numerical  method  of  Chen  and  Chen  (1982,  1984),  devel¬ 

oped  initially  for  elliptic  equations,  is  then  revised  and  ex¬ 
tended  to  solve  the  five  transport  equations  (IV-1)  for  4>  = 
U , V , W , k  and  e  with  a  guessed  pressure  field.  This  is  described 
in  Section  IV. 2.  The  continuity  equation  is  used  to  obtain  equa¬ 
tions  which  enable  the  determination  of  the  pressure  and  pres¬ 
sure-correction  fields  from  the  velocity  field.  This  procedure 
is  described  in  Section  IV. 3.  The  treatment  of  the  boundary 
conditions  is  described  in  Section  IV. 4,  and  the  overall 
numerical  solution  algorithm  is  then  outlined  in  Section  IV. 5. 

IV. 1.  Grid  Generation 

As  described  in  Section  1 1 1 . 1 ,  a  general  curvilinear  co¬ 
ordinate  system  can  be  constructed  by  numerically  solving  equa¬ 
tions  (  1 1 1-2 1  )  ,  i.e., 
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h  h  h  1  '  h  ' 
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V2X2  =  _ <L_  (  V|3v 
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( IV-3) 


V2x3  =  — 1 _ 

h,h„h^  .  3  '  h.,  ’ 

1  2  3  3x  3 
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( IV-3a ) 


are  the  metric  coefficients  in  the  chosen  orthogonal  coordi¬ 
nates  x1  and  g1^  is  the  conjugate  metric  tensor  in  the  trans¬ 
formed  coordinates  (  £  ,  n  ,  ?  ) .  The  control  functions  f1  can  be 
chosen  to  generate  the  desired  grid  concentration.  Note  that  the 
orthogonal  coordinates  x1  in  the  physical  plane  can  be  selected 
arbitrarily.  The  equations  for  two  cases  are  given  below. 

(i)  Cartesian  Coordinates  in  the  Physical  Plane 

When  Cartesian  coordinates  are  used  to  describe  the  three- 
dimensional  computational  domain,  a  body-fitted  grid  system  can 

1  9  O 

be  generated  from  equation  (IV-3)  with  (x  ,  x  ,  x  )  =  (x,y,z)  and 
(h1 ,h2,h3)  =  (1,1,1),  i.e.  , 
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( IV- 5  ) 


the 


Hor  a  two-dimensional  flow  field  in  the  xy-plane, 
rning  equations  (IV-4)  can  be  simplified  by  letting  s  = 


z ,  so 
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( IV- 6  ) 


(  IV-  7 ) 


(ii)  Cylindrical  Polar  Coordinates  in  the  Physical  Plane 


For  applications  related  to  axisymmetric  bodies  and  ship 
forms,  it  is  more  convenient  to  use  the  cylindrical  polar  coordi¬ 
nates  instead  of  the  Cartesian  coordinates.  The  governing  equa- 
tions  for  this  transformation  are  obtained  by  letting  (x  ,  x  , 


x3)  =  (x,r,0)  and  (h1,h2,h3)  =  (1,1, r)  in  equation  (IV-3),  i.e., 
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(iii)  Choice  of  Coordinates  and  Numerical  Solutions 

The  above  equations  can  be  solved  numerically  in  the  trans¬ 
formed  domain  £ 1  when  appropriate  grid-control  functions  f1  are 
selected  and  proper  boundary  conditions  are  specified  on  all 
boundary  surfaces.  Depending  upon  the  particular  application, 
equations  (IV-4)  or  (IV-8)  have  been  used  previously  to  generate 
various  body-fitted  coordinates.  A  review  of  these  algorithms 
has  been  made  recently  by  Thompson,  Warsi  and  Mastin  (1982).  The 
present  formulation  provides  greater  flexibility  and  generality 
in  the  selection  of  the  coordinate  systems  in  the  physical  plane. 


'14 


Insofar  as  the  problem  of  grid  generation  is  concerned,  the 
choice  of  the  basic  orthogonal  coordinate  system  used  to  describe 
the  boundaries  of  the  computation  domain  in  the  physical  plane  is 
a  matter  of  convenience.  For  example,  the  flow  on  an  axisym- 
metric  body  is  best  described  in  cylindrical  polar  coordinates 
although  a  three-dimensional  formulation  with  Cartesian  coordi¬ 
nates  can  also  be  employed.  If  the  boundaries  in  the  physical 
plane  are  all  straight,  it  may  be  more  convenient  to  start  with 
the  Cartesian  coordinates.  For  three-dimensional  bodies,  such  as 
ship  hulls,  both  the  Cartesian  and  the  cylindrical-polar  coordi¬ 
nates  can  be  used.  If  the  equations  of  motion  are  transformed 
completely  so  that  the  velocity  components  are  in  the  direction 
of  the  transformed  coordinates,  the  choice  of  the  orthogonal 
coordinates  used  in  the  physical  plane  is  immaterial.  However, 

I  if  only  a  partial  transformation  is  used,  with  the  velocity  com¬ 

ponents  left  in  the  basic  orthogonal  system,  as  in  the  present 
formulation,  then  the  choice  of  the  basic  coordinate  system  be¬ 
comes  important  and  is  dictated  by  the  numerical  algorithm  used 
j  to  solve  the  flow  equations.  The  problem  is  illustrated  in  Fig. 

3,  which  shows  a  numerically-generated  grid  for  a  ship  section. 
It  is  clear  that  the  use  of  Cartesian  components  which  do  not 
rotate  with  the  n ,  c  coordinates  will  lead  to  serious,  and  perhaps 
|  insurmountable,  numerical  difficulties.  On  the  other  hand,  since 

the  velocity  components  in  polar  coordinates  rotate  with  the  n , t 
directions,  the  numerical  problems  are  avoided.  Although,  as 
noted  above,  the  need  for  an  a  priori  choice  of  the  basic 
I  coordinates  in  the  physical  plane  for  each  class  of  problems  can 

be  avoided  by  a  complete  transformation  of  the  equations  into  the 
numerical  cocord i nates ,  the  additional  computer  storage  required 
for  practical  three-dimensional  flows  appears  to  be  prohibitive 
I  at  the  present  time,  and  the  authors  are  not  aware  of  successful 

solutions  using  such  coordinates.  However,  the  partial  trans¬ 
formations  used  here  enable  solutions  to  be  obtained  in  a  cost- 
effective  manner  with  modest-size  computers. 
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For  applications  to  axisymmetric  bodies  and  elongated  three- 
dimensional  bodies  with  one  or  more  planes  of  symmetry,  such  as 
airplane  fuselages  and  ship  forms,  it  is  convenient  to  construct 
coordinates  such  that  axial  cross-sections  of  the  body  form  one 
set  of  coordinate  surfaces,  i.e.  to  select  t(x),  where  x  is 
measured  alone  the.'  body  axis.  This  results  in  some  simplifica¬ 
tion  of  the  equations  of  motion  and  also  facilitates  the  specifi¬ 
cation  of  the  computation  domain  and  boundary  conditions.  Other 
advan tages  include  the  possibility  of  constructing  orthogonal  or 
near  1 y-orthogonal  coordinates  in  the  cross-sectional  planes,  a 
continuous  evolution  of  the  coordinates  from  the  tail  of  the  body 
into  the  wake,  an  i  the  ease  of  interpretation  of  the  results.  An 
obvious  disadvantage  of  us  i  ••  generalized  coordinates,  which  is 
>f ten  overlooked ,  is  that  cons i durable  interpolation  in  the  com¬ 
muted  results  may  be  required  t  •  o. der  to  make  comparisons  with 
experimental  data  obtained  in  a  convenient,  and  usually  ortho- 
jona L ,  laboratory  coordinate  system. 

;  he  present  coord  nets  aeration  method  is  illustrated  by 

means  of  two  examples,  nas-ly  an  axisymmetric  body  and  an  elon¬ 
gated  ship-like  body  win  cross-sect  ions  are  everywhere  elliptic 
with  an  axis  ratio  re  1:1.  rn  both  cases  the  flow  region  of 

interest  is  the  thick  boundary  layer  over  the  stern  and  its  evol¬ 
ution  into  the  wake. 

To  genera'  •>  t  he  coordinates  for  an  axisymmetric  body  (After- 
Pod/  S  of  Muang  ->t  al  .  ,  I9«fll  ,  the  constant-^  stations  are  chosen 
t.o  be  the  const  ant-x  pianos  (i.e.,  the  transverse  sections).  The 
r  i  r  at.  .t.  di  e,  =  1,  is  located  at  X  =  0.6L  while  the  last 

station,  ?  ~  60,  is  located  in  the  far  wake  at  X  =  13.20L,  where 

r,  is  the  body  length  and  x  is  the  axial  distance  from  the  nose. 
Tn  the  radial  direction,  the  outer  boundary  is  cylindrical  and  is 
placed  at  W  -  n.  ti  171,,  which  corresponds  to  n  =  19  in  the 
transformed  domain.  tinder  this  arrangement,  equations  (IV-9) 
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The  control  function  f-'-  is  therefore  determined  by  the  desired 
distribution  of  the  axial  stations.  With  f^  specified,  equation 
(IV-11)  yields  the  distribution  of  points  in  the  radial 
direction,  r(5,n).  To  obtain  the  desired  grid  distribution  in 
this  direction,  the  control  function  f*  is  prescribed  by 
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where  C  =  15,  5  =  42  ,  and  F.  and  Fn  are  determined  by  the  node 
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distributions  at  the  initial  and  final  stations  as 
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(  IV- 13b) 
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and  Fc  is  given  by  a  linear  combination  of  FA  and  F^: 


(  IV- 13c) 
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Note  that  equation  (IV-11),  together  with  the  f  defined  in 
equation  (IV-13),  reduces  to  the  following  homogeneous  equation 
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which  is  equivalent  to  the  two-dimensional  formulation  in  Carte- 
sian  coordinates  (x,r)  with  control  functions  f  and  F  .  There¬ 
fore,  the  same  grid  distributions  can  be  generated  in  both  the 
cylindrical  and.  Cartesian  coordinates  if  the  control  function  f 
is  replaced  by  Fz  in  the  latter  formulation. 

Equation  (IV-14),  with  the  control  functions  specified  by 
equations  (IV-12)  and  (IV-13),  can  be  solved  by  any  convenient 
numerical  method.  Here,  we  have  used  two  different  finite- 
difference  techniques,  namely,  a  central-difference  scheme  and  an 
exponential  scheme. 

In  the  central-difference  scheme ,  the  discretized  form  of 
equation  (IV-14)  can  be  written 


r  =  { (g1 *+  0 . 5  f  ^ )  ,  _r 
s  f  n 


f  _ - c . ,  +  ( g 1 1  -  O.sf1),  r  , 

5  i  f.+  l  ,  n  £  ,  n  £-1  ,  n 


+  °-5F\,n^,n.i+  ^22"  °‘5F2)C,n^,n-l 

+  °’59£,n  (re  +  l,*iH+  r£-l,n-l~  r£-l  ,  n  +  l~r£  +  l  ,  n-1 } } 
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/12(q  +  g  )^n 


( IV- 1 5 ) 


with  2<£<59,  2<n<18.  With  the  Neumann  conditions  specified  at  the 
boundaries,  this  system  of  algebraic  equations  was  solved  by  an 
SOR  algorithm  with  an  over- re laxa t ion  factor  of  1.8.  The  results 
are  shown  in  Fig.  4  by  the  dotted  lines.  The  bounda  ry- f  i  t  ted 
coordinates  thus  obtained  evolve  smoothly  from  the  body  into  the 
wake,  and  give  the  desired  grid  distribution  in  the  thick 
boundary  layer  region  over  the  tail  of  the  body. 

•  1  2 

We  note  that,  with  nonzero  f  and  F  ,  equation  (IV-14)  is, 

in  fact,  a  convective-transport  equation  with  convective 

111  222 

velocities  or  cell  Reynolds  numbers  -  f  /g  and  -  F  /g  "  in  the 


4  8 


transformed  plane.  It  is  well  known  that  a  negative  F  (C,n) 
value  (i.e.,  positive  convective  velocity)  will  result  in  an 
upwind  bias  from  smaller  constant-n  lines,  and  therefore  pull  the 
grid  lines  toward  the  body  surface.  On  the  other  hand,  a 


positive  value  results  in  the  expansion  of  grid  lines  near  the 


body  surface.  For  large  values  of  the  cell  Reynolds  numbers 


i .  e  .  , 
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f  /g  or  F  /g  >2),  the  central-difference  scheme 


becomes  inaccurate  and  unstable  due  to  the  negative  coefficients 


(i.e.,  (g 1 1  +  0.5fl)  <  0  or  { g ^  +  0.5F^)  <  0)  appearing  in  the 


discretized  equation  (IV-15). 


In  order  to  avoid  the  instability  of  the  central-difference 
scheme  and  improve  the  accuracy,  the  exponential  scheme  of 
Spalding  (1972)  was  also  employed  to  provide  an  automatic  upwind 
shift.  With  this  scheme,  the  governing  equation  (IV-14)  can  be 


written  as 


11  22 

( 2g  a  coth  a  +  2g  b  coth  b)r  rr 

s  t  n  s  /  n 


=  (2gUa  csch  %  +  Un+  earc_^n) 


,_22,  ,  .  .  ,  -b  ,  b  . 

f  <2g  b  csch  b)^n(e  r^n+1+  e  r^^) 


+  2gS,n  (r£+l,n+l+  rC-l,n-l“  r5-l,n+l  r^+l,n-l) 


( IV- 1 6 ) 


with 


2a  =  - 


( IV-16a) 


5=constant 


2b  =  - 


( IV- 16b) 


The  above  system  of  algebraic  equations  was  solved  by  the  tri- 
<1  iayonal-matr  ix  algorithm  with  an  over-relaxation  factor  of 
1.8.  The  coordinates  thus  obtained  are  shown  by  the  solid  lines 
in  Fig.  4 . 

It  can  be  seen  that  the  coordinate  lines  obtained  by  both 

schemes  agree  quite  well  everywhere.  However,  the  discrepancy 

2 

will  increase  when  larger  values  of  F  (£,n)  are  used.  Note  that 
the  central-difference  scheme  yields  consistently  lower  values 
due  to  the  stronger  upwind  influence  from  smaller  constant-n 
lines.  This  eventually  leads  to  the  occurrence  of  a  minimum  in 
the  interior  of  the  domain  as  the  values  of  F  become  large.  On 
the  other  hand,  the  exponential  scheme  assures  a  one-to-one 
mapping  between  the  physical  and  transformed  planes  since  all 
influence  coefficients  in  equation  iIV-16)  are  positive. 

In  order  to  demonstrate  the  grid-generation  technique  for  a 
three-dimensional  body,  we  consider  the  3:1  elliptic  body  of 
Huang  et  al.  (  1983).  In  this  example,  the  outer  boundary  is  a 
cylindrical  surface  of  radius  R  =  0.8137L,  which  corresponds  to 
n  =  19  in  the  transformed  plane.  The  constant-^  stations  are 

again  chosen  to  be  the  constant-x  planes,  with  the  initial  (£  =  1) 
and  final  (£=60)  stations  located  at  X  =  0.5L  and  X  =  16.25L, 

respectively.  In  the  circumferential  direction,  seven  stations 
are  used,  with  C=2  and  6  corresponding  to  8  =  0°  and  90°,  respec¬ 
tively.  The  symmetry  conditions  at  9  =  0°  and  90°  are,  respec¬ 
tively. 


r  (  £  ,  n  ,  1 ) 

=  r  (  £  ,  n  ,  3  )  , 

9  (  £  ,  n  ,  1 ) 

=  -  9  (  £ 

f  n  ,  3  ) 

r  (  £  ,  n  ,  7 ) 

=  r  (  £  ,  n  ,  5  )  , 

w  (  £  ,  n  ,  7 ) 

=  Ti-e  ( £ 

,n,5) 

at  9  =  0 
at  6  =  1 


(  IV- 1 7 ) 


Under  this  arrangement,  equations  (1V-8)  become 


11  22  33  12  13  „  23 

g  r  +  g  r  +  g  r  +  2g  r  +  2g  r  +  2g  r 


12  3  1 

+  f  r5+  f  rn  +  f  r^  =  ~ 


uv-y) 


11Q  22  -  3  3  .  „  1 2  Q  ^  13.  ^  2  3  Q 

g  0CC+  g  %+  g  2g  8^+  2g  e?;+  2g  0n4 

12  3 

+  f  6+  fz0  +  fJ0  =  0 
5  n  5 

The  control  functions  f*-  and  defined  by  equations  (  I V—  1 2  >  and 
(IV-13)  are  again  employed  here,  while  f^(f)  is  determined  by  the 
circumferential  grid  distribution  on  the  outer  boundary  as 


,3,  .  33  CC 

f  ( e)  =  -  g  g — 

5  5  =  1  ,  n  =  1 9 


( IV- 1 9 ) 


In  order  to  facilitate  the  use  of  wall  functions  in  the 
turbulence  model,  Neumann  boundary  conditions  are  specified  on 
all  boundaries;  i.e.,  n  =  1,19;  5  =  2,6;  of  the  contant-5  sta¬ 
tions  to  ensure  local  orthogonality  of  these  coordinates.  Equa¬ 
tions  ( IV- 1 8 )  are  then  solved  by  a  finite-difference  method. 

Because  the  central-difference  scheme  used  in  the  previous 

1  2 

example  indicated  a  lack  of  convergence  for  large  f  ,  f  and/or 
-3 

f  ,  only  the  exponential  scheme  described  above  was  employed 
here.  After  some  manipulation,  the  discretized  equations  become 


(2g^a  coth  a  +  2g22b  coth  b  +  2g^c  coth  c) 


5  ,  t ,  C  5  ,  n  ,  C 


-  ( 2cj  a  csch  a)  (e  ^ d>  +  ) 

y  d  6  ,  n  ,  £ ' e  ?5+l,n,f  e  *6-l,n,C 

i  (2g22b  csch  b)r  r(e  h<J>r  .  +  eh4>  -  ,  ,) 

5  ,  n  ,  t  5 , n+ 1  ,  t  5  ,  n- 1  ,  f 

3 1  — .  £•  c 

+  (2g  c  csch  c)  ,  (e  r  „  o  >}> .  r  .) 

5  ,  n  ,  5  6  ,  n  ,  5  +  1  6  ,  n  ,  ;  -  1 


+  2a  (A  +  a  -4) 

J  4  ,n  ,  r/  *4+1  ,  n  +  1  ,  r,  ,  n-1  ,  r,  *4+1  ,  n-1  ,4 


*5-1,  n+1,51  +  2q  5,n,c  U5  +  l  ,n,e+l+  *5-l,n,C-l 
*4+l,n,c-l_  *5-1, n  ,5  +  1*  +  2g5  ,  n  ,  e  (  *4  ,  n  +  1  ,  5  +  1 


4*  til  —  <fi  —  d>  ) 

v  5  ,  n  - 1  ,5-1  5,  n+1,  c-1  5,n-l,5  +  l 


( IV-20 ) 


where  4>  represents  either  r  or  9,  and 


2a  =  - 


11, 

g  5  ,  n  ,  c 


‘5  15  =  constant. 


( IV- 20 a  ) 
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2  b  = 


FA(n>  5<5a 

Fc(f'r)  V^b 


Fp(n) 


4  >  5 

B 


( IV-20b) 
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g  t. ,  n  , 


for  $ 


2c  =  - 


4  =  1  ,  n  -  1  y 


IV-20c) 


P  (  n )  = 

A  r 


n  !  5  =  1  ,  5  =  6 


tM(n)  =  r-~ 


4  =  1 0  ,  '  -  h 


t(,  U,n) 


4  -  5 )  P  +  (4-4  )  F  ]  /  (  4  -  4  ) 
It  A  A  H  B  A 


Not*1  that  t  he  coil  r  i  f..jru  ’  i  on;'  F"  and  1  *  can  be,  in  general, 
i  i.i  n  -t  ions  of  ■  ,  n  and  4,  although  simpler  expressions  are  used 
her'.'.  Fat  more  complicated  geometries,  the  control  functions  may 
!»•  determine';  by  the  node  distributions  on  all  boundary  surfaces 


and  an  interpolation  between  them.  Detailed  discussions 

regarding  the  selection  of  control  functions  can  be  found  in, 
among  others,  Thompson  (  1982)  and  Thompson,  Warsi  and  -v  a  s  t  '  n 
(  1982  )  . 

With  the  Neumann  boundary  conditions  specified  on  all 
boundaries  of  the  constant-^  stations,  equations  ( IV-20 )  were 
solved  by  the  tridiagonal-matrix  algorithm  with  an  over¬ 
relaxation  factor  of  1.4.  The  coordinates  thus  obtained  are 

nearly  orthogonal  at  each  cross-section,  and  also  provide  the 

desired  expansion  in  the  thick  boundary  layer  region  over  the 

tail.  Figure  5  shows  several  views  of  the  coordinates  in  the 
stern  region. 

In  the  present  study,  the  control  functions  f1  were  either 
specified  beforehand  or  evaluated  numerically  from  the  prescribed 
distributions  of  boundary  nodes.  In  most  applications,  the 
geometric  coefficients  g  ^  j  ,  g1^  and  b;?  were  evaluated  using 

central-difference  formulae.  However,  more  accurate  exponential 
expressions,  derived  from  the  local  analytic  solutions,  can  be 
employed  to  evaluate  these  geometric  coefficients  when  the  grid 
expansion  ratios  become  large.  Note  that  coordinate  systems 
generated  by  central-difference  and  exponential  schemes  will,  in 
general,  give  different  geometric  coefficients  even  if  the  same 
control  functions  are  employed. 

Since  the  geometric  coefficients  and  the  control  functions 
appear  in  the  continuity  and  the  convective-transport  equations, 
numerical  errors  arising  from  the  grid-generation  technique  are 
not  entirely  decoupled  from  errors  in  the  flow  calculations.  In 
order  to  reduce  this  error  propagation,  the  coordinates  generated 
by  the  exponential  scheme  were  adopted  in  all  calculations 
presented  here. 

If  other  grid-generation  techniques,  such  as  conformal 
transformations,  are  employed  to  generate  the  coordinates,  it  is 
sti  11  possible  to  use  the  convoc t i ve- 1 ranspor t  equations  (  TV-1  ) 


derived  earlier,  provided  the  control  functions  f1  are  evaluated 
numerically  from  equation  (A-57),  i.e., 

f1  =  4  ~T  )  (IV-21) 

d  (,  1 

In  this  fashion,  numerical  errors  in  the  flow  calculations  can  be 
decoupled  from  those  introduced  by  the  grid-generation 
technique.  However,  the  use  of  equation  (A-57)  for  f1,  which 
involves  second  derivatives  of  the  coordinates,  may  result  in 
significant  numerical  errors  in  the  fluid-flow  calculations. 
Therefore,  the  control  functions  used  here  were  read  directly 
from  the  grid-generation  routines,  rather  than  evaluating  them 
numerically  from  equation  (A-57). 

IV. 2  Finite-Analytic  Method  for  the  Transport  Equations 

Because  of  the  absence  of  the  term  <t>  in  the  convective 

S»  ^ 

transport  equations  (IV-1),  these  equations  are  parabolic  in 
the  f,~  direction  but  remain  elliptic  in  the  nt  plane.  In  order 
to  correctly  predict  the  elliptic  nature  of  the  flow  at  each 
cross-section,  the  f i n i t :-analy t ic  method  of  Chen  and  Chen  (1982, 
1984),  has  been  revised  and  extended  to  solve  the  present 
part ial ly-parabol ic  transport  equations.  The  solution  procedure 
is  described  in  the  following. 

In  the  finite-analytic  approach,  equations  (IV-1)  in  each 
rectangular  numerical  element,  A£>An=Ac=l,  are  locally  linearized 

as 


2  2  3  3 

’  p  n  n  J  p 
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+  2  (  B  ,  ) 
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(  IV-22) 


w  i  t.  h 


IVii  u  _  ^n-1 

r  '  r  D  M) 


+  (S.)D 


( r> , )  ,,{p 


(  I V-2  2a ) 


where  the  subscript  P  denotes  the  center  node  p  of  the  element 
shown  in  Fig.  6  and  t  is  the  time  increment.  Thus,  the  time-  and 
derivatives  are  approximated  by  backward-difference  formulae, 
subscript  (J  denoting  the  upstream  nodal  value,  and  superscript 
(  n— 1 )  denoting  the  value  at  the  previous  time  step.  If  we  intro¬ 
duce  the  coord i nate-st re tch i ng  functions 


★ 

n 


( IV- 2 3 ) 


equation  (IV-22)  reduces  to  the  standard  two-dimensional  con¬ 
vective-transport  equation  described  in  Chen  &  Chen  (1982,  1984), 

i  .  e  .  , 


1  *  *+  ♦  *  2A4>  *+  2B  $  *+  g 

C  C  n  n  c  n 


( IV-24 ) 


with 
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( IV-24a ) 


for  a  local  element  with  dimensions 
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( IV-24b ) 


By  specifying  a  combination  of  exponential  and  linear 

boundary  functions,  which  are  derived  from  the  natural  solutions 

★ 

of  the  governing  equations,  on  all  four  boundaries  n  =  ±k 

★ 

and  C  =  1  h  of  each  local  element,  i.e.. 


(  IV-25a ) 
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(  IV- 25b ) 

( IV-25C) 

(  I V-2  5d  ) 


whore  a,  h  and.  >:  arc-  constants,  equation  (IV-24)  can  be  solved 
analytically  tor  each  element  by  the  method  of  separation  of 
variables  or  any  other  analytic  technique.  Details  of  the 
solution  procedures  are  described  in  Chen  6.  Chen  (1982  ,  1984). 

When  the  local  analytic  solution  thus  derived  is  evaluated  at  the 
central  node  P  of  the  element,  the  following  nine-point  finite- 
analytic  algebraic  equation  is  obtained. 


h  ~  c  t  +  c  +  C  +  c  $  +  C  $  +  C  ^ 

9P  Nr  NK  NWfNW  "Kir  Si:  SVT  SW  ECTEC  WCVWC 
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( IV-26 ) 
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summation  in  E2  can  be  avoided  by  considering  the  asymptotic 
expressions  of  PA  and  PR  based  on  the  theory  of  characteristics, 
i  .  e  .  , 


|  Ak  j  >  | Bh |  :  PA  =  0,  PR=  1  -  j  Bh/Ak | 
|  Ak |  <  | Bh |  :  PR  =  0  ,  Pft=  1  -  |Ak/Bh| 


(  IV- 28 ) 


Since  the  downstream  influence  is  negligible  at  large  cell  Rey¬ 
nolds  numbers,  the  above  approximations  do  not  introduce  a 
significant  error  in  the  solution,  but  the  computing  time  is 
greatly  reduced. 

By  substituting  the  nonhomoyeneous  term  g  from  eguation  (IV- 
22a)  into  equation  (IV-26),  an  eleven-point  finite-analytic 
formula  for  unsteady,  three-dimensional,  partially-parabolic 
equations  can  be  obtained  in  the  form 

^P  =  nrr~  *cne*ne+  cnwi*>nw+'  csei<>se+  CSVI*SVI 
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where  the  subscript  nb  denotes  neighboring  nodes  ( NE : northeast , 
NW: northwest ,  etc). 
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It  is  seen  that  $  depends  upon  all  the  eight  neighboring 
nodal  values  in  the  n-c  plane  as  well  as  the  values  at  the 

upstream  node  U  and  the  value  at  the  previous  time  step  (n-1). 
Since  the  finite-analytic  coefficients  are  all  positive,  the 

ellipticity  of  the  flow  in  each  cross-section  is  represented 
accurately.  Note  that  is  not  directly  dependent  upon  the 
downstream  nodal  value  This  parabolic  feature  enables  us  to 

employ  a  marching  procedure  along  the  £-  direction  as  well  as  in 
time.  However,  the  flow  is  not  truly  parabolic  in  the  £- 

direction  since  the  longitudinal  pressure  gradient  P  ,  which 
appears  in  the  source  functions,  introduces  the  downstream 
influence  of  the  elliptic  pressure  field.  Due  to  this  indirect 
pressure  transmission,  several  sweeps  in  the  £-  direction  are 
needed  to  obtain  a  fully-converged  pressure  and  velocity  field. 

Since  equations  (  I V—  29)  are  implicit,  both  in  space  and 
time,  at  the  current  station  of  calculation,  their  assembly  for 
all  elements  results  in  a  set  of  simultaneous  algebraic 

equations.  These  equations  are  solved  by  the  tridiagonal-matrix 
algorithm.  Since  it  is  not  necessary  to  obtain  a  fully-converged 
intermediate  solution  for  steady  flows,  only  six  line-by-line 
internal  iterations  are  used  during  each  global  sweep. 
Furthermore,  the  f i n i te-a na 1 y t i c  coefficients  appearing  in 
equations  (IV-29)  are  not  updated  during  these  internal 
iterations  for  economy  of  computation  time. 

IV. 3  Solution  of  the  Continuity  Equation,  and  Pressure  and 

Velocity  Corrections 

If  the  pressure  is  known  a  priori,  equations  (IV-29)  can  be 
employed  directly  to  solve  the  five  part ial ly-parahol ic  convec¬ 
tive-diffusion  equations  ( I V— 1 )  for  U,V,W,k  and  e  .  However,  in 
most  practical  applications,  the  pressure  is  not  known  and  must 
be  determined  by  requiring  the  velocity  field  to  satisfy  the 
equation  of  continuity.  Since  a  direct  method  for  the 

simultaneous  solution  of  all  six  equations  is  not  feasible  with 
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present  computer  capacity,  it  is  necessary  to  convert  the  equa¬ 
tion  of  continuity  into  an  algorithm  for  the  calculation  of  the 
pressure  field  appearing  in  the  momentum  equations.  In  this 
study,  the  SIMPLER  algorithm  of  Patankar  (1980)  has  been  modified 
and  extended  to  update  the  pressure  field.  Details  of  this  al¬ 
gorithm  are  outlined  in  the  following. 

In  the  prt sent  algorithm,  the  staggered-gr id  system  intro¬ 
duced  by  Harlow  and  Welch  (1965)  in  their  MAC  method  is 
adopted.  Figure  6(b)  shows  the  locations  of  the  nodes  for  U,V,W, 
and  P  in  the  staggered  grid.  The  turbulence  quantities  k  and  e 
are  evaluated  at  the  pressure  nodes.  The  dashed  lines  represent 
the  control  volume  faces,  and  the  pressure  is  calculated  at  the 

center  of  the  control  volume.  For  convenience  IK,  V  ,  W  and 

d  n  e  *^p 

in  Fig.  6(b)  are  assigned  the  same  index,  i.e.,  they  are  denoted 

by  U£n£'  VCnt'  W£n;  and  p£nt'  resPect ively .  It  should  be 

recalled  here  that  the  velocity  components  U,V,  and  W  remain  in 
the  longitudinal,  radial  and  circumferential  directions.  In 
other  words,  they  are,  in  general,  neither  perpendicular  to  the 
control  surfaces  nor  in  the  direction  of  the  coordinate  lines. 
However,  due  to  the  choice  of  the  cylindrical  polar  coordinates 
in  the  physical  plane,  the  velocity  components  do  not  become 
parallel  to  the  control  surfaces. 

In  the  staggered  grid,  the  eleven-point  finite-analytic 
formulae  (equation  (IV-29))  for  the  momentum  equations  yield  the 
velocity  components 
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where  ,  Cn  and  Ce  are  the  f  i  n  i  te-analy  t  ic  coefficients  Cp 
evaluated  at  the  staggered  velocity  nodes  d,  n  and  e  in  Fig.  6. 
Note  that  the  above  equations  contain  the  pressure-gradient  terms 
inside  the  source  functions.  An  equation  for  this  unknown 
pressure  field  may  be  obtained  as  follows. 


If  we  decompose  the  actual  velocity  field  (UfV,W)  in  the 
momentum  equations  (IV-30)  into  a  pseudovelocity  field 
(U,  V,  W)  plus  the  pressure  gradient  terms  contained  in  the 
source  functions,  i.e.. 
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so  that  the  pseudovelocities  contain  no  pressure  terms,  then  an 
equation  for  pressure  can  be  derived  by  requiring  the  velocity 
field  to  satisfy  the  discretized  equation  of  continuity  (IV-2), 
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( IV- 32) 


The  resulting  pressure  equation  will  contain  many  pressure  nodes 
(see  Muraoka  (1979),  for  example)  if  nonorthogonal  coordinates 
are  employed.  It  is  therefore  desirable  to  simplify  the 
algorithm  by  introducing  modified  pseudovelocities  (U,  V,  W)  by 
decomposing  the  velocity  components  as  follows: 
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(  IV-33a ) 


The  modified  pseudovelocities  (u,  V,  W)  defined  above  still 
contain  part  of  the  pressure-grad ient  terms  if  the  coordinate 


system  is  nonorthogonal  (i.e.,  b|  *  0  for  i  *■  j  ,  see  equation 
IV-31).  These  pressure-gradient  terms  can  be  evaluated  from  the 
pressure  field  known  at  the  previous  time  step  or  iteration  with¬ 
out  losing  any  accuracy  or  generality.  If  we  require  the 
velocity  field  to  satisfy  the  equation  of  continuity  (IV-32),  a 
simpler  pressure  equation  can  be  derived  in  terms  of  the  modified 
pseudovelocities  (U,  v,  W)  .  Note  that  eighteen  velocity  compo¬ 
nents  are  involved  in  the  equation  of  continuity  (IV-32)  for  each 
small  control  volume.  However,  due  to  the  staggered  grid  system 
employed  here,  only  six  of  these,  namely,  ,  Uu ,  Vn ,  Vg ,  We ,  and 
Ww,  can  be  obtained  directly  from  the  governing  equations  (IV- 
30).  It  is,  therefore,  necessary  to  approximate  the  remaining 
twelve  by  interpolations.  A  simple  linear  interpolation  is  used 
here  to  evaluate  these  values  from  the  velocity  field  known  at 
the  previous  time  step  or  iteration,  so  that  the  continuity  equa¬ 
tion  becomes 

(bju).  -  ( bj  U)  +  ( b^V)  -  ( b? V)  +  (b^W)  -  (b*W)+  D  =  0  (IV-34) 

Id  lu  2n  2s  3  e  3wl 

where 

D,  =  (biv  +  bjw)  .  -  (b^v  +  bjw)  +  (b?u  +  b,w)  -  (b?u  +  b^w) 

1  2  3d  2  Jul  ini  3s 

+  (b^U  +  b2V)e  ~  (blU  +  b2V)w  ( IV- 34a ) 

is  the  mass  source  obtained  from  the  velocity  field  at  the 
previous  time  step  or  sweep.  An  equation  for  pressure  is  then 
derived  by  substituting  equations  (IV-33)  into  (IV-34),  i.e., 

apPP  =  adPD  +  aupU  +  anpNC  +  aspSC  +  aepEC  +  awpWC  ”  D  (IV-35) 
where 

ad  =  (  bi  dd 
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(  IV-35a ) 


The  modified  pseudovelocities  (Uf  V,  W)  contain  the  neigh¬ 
boring  nodal  values  of  velocity,  source  functions,  and  part  of 
the  pressure-gradient  terms.  All  of  them  can  be  evaluated  from 
the  information  known  at  the  previous  time  step  or  iteration. 
Therefore,  apart  from  the  interpolations  for  ,  the  pressure 

equation  (1V-35)  is  still  an  exact  algebraic  representation  of 

the  equation  of  continuity  (IV-2).  In  this  fashion,  the  pressure 
field  can  be  obtained  directly  from  an  estimated  velocity  field, 
and  the  slow  convergence  usually  encountered  in  the  SIMPLE 
algorithm  is  avoided. 

Although  the  queue. -i  pressure  field  can  be  updated  directly 
by  equation  (  1 V—  3  S i  ,  ;  p  practical  applications  the  new  pressure 

field  may  produce  i  i  i  t  y  t  1  e  1  d  which  does  not  satisfy  the 
equation  of  con  *-  :  n  .  .  ‘  .  .  »'•  iterative  procedure  is  therefore 

required  to  cor  n  .  ■  t  t  r  i  :  •  ■  r  r  .  is  velocity  field  to  achieve  more 

rapid  convergence.  it,  thin  Mudy,  a  velocity-correction  formula, 


similar  to  that  used  in  the  SIMPLE  algorithm,  is  derived  in  terms 
of  the  pressure-corrections. 

If  we  denote  the  imperfect  velocity  field  resulting  from  an 
imperfect  pressure  field  p*  by  (U*,V*,  W*  )  ,  then  the  discretized 
momentum  equations  (IV-33)  can  be  written  as 

*  ~  *  *  * 

Ud  =  Ud  "  de{pD'  PP> 

Vn  =  Vn*“  dn(pNC-  PP^ }  (IV-36) 

★  A  *  ★  ★ 

We  *  We  -  de(PEc‘  Pp> 


In  order  to  improve  the  guessed  pressure  field,  such  that  (U*, 
V* ,  W* )  will  eventually  satisfy  the  equation  of  continuity,  one 
needs  to  know  how  the  velocity  components  respond  to  a  change  in 
the  pressure  field.  Such  a  relation  can  be  obtained  by  sub¬ 
tracting  equation  (IV-36)  from  equation  (IV-33),  i.e., 
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where  p' 

1  =  p  - 

P* 

is  the 

pressure  correction, 

and  U-U*,  V-V* 

and 

W-W*  are  the  corresponding  velocity  corrections.  If  we  require 
the  velocity  field  to  satisfy  the  equation  of  continuity  (IV-34), 
an  equation  for  the  pressure  correction  p'  can  be  derived. 
However,  due  to  the  implicit  nature  of  the  velocity  corrections 
arising  from  the  pseudovelocities,  the  resulting  pressure 
correction  equation  would  involve  the  pressure-corrections  at  all 
grid  points.  It  is  not  necessary  to  retain  such  a  complicated 
formulation  because  both  the  pressure-  and  velocity-correct  ions 
are  zero  in  the  final  converged  solution.  Since  both  the 


pressure- 


become  trivial  when  the 


and  velocity- 

solution  converges,  it.  is  possible  to  omit  that  part  of  the 

'  -  *  ~  *  *  *  -  * 

velocity-corrections,  U  -  U  ,  V  -  V  and  W  -  W  ,  which  repre¬ 

sents  the  indite:*-  influence  of  velocity-corrections.  With  this 
approximation,  the  ve  loci  ty-eorreet.  ions  are  expressed  explicitly 
in  terms  of  tin  p  r  esse  i  >  — ■  r.  ir  re.- 1  1  .vis  cs 
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By  requiring  the  velocity  field  'o  satisfy  the  equation  of  con¬ 
tinuity  (  IV-30  )  ,  a  pressur.i-co:  ;  oc  *- i  on  equation  is  then  obtained 
in  the  form 


appp=  adpD  +  aupllf  anp'-:f  :;sp,.r+  depKC+  awpWC~  D 
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when*  the  coefficients  apf  a^  ,  etc.  are  as  defined  in  equation 
MV-3Sa).  Note  that  the  pressure-correction  equation  is  similar 
to  the  pressure  equation  (IV-35).  Although,  unlike  the  pressure 
e  'nation,  the  pressure-correction  equation  is  not  exact,  the 
ippro.x  imat  l  ons  made  will  influence  only  the  rate  of  convergence 
but  not  t  t!<>  i  i  no  I  converged  solution. 

It  s  1 1  >  i !  t  be  t'ema  rkeo  here  that  both  the  pressure  and  p  res¬ 
in  t  —  - 1  ,  e(  -  ■  i .  ■  n  egu  a*  ions  t  r  c  elliptic.  In  other  word  s  ,  the 

downst  ream  pressure  p()  and  pressure-correct.  ion  will  affect  the 


pressure,  and  thus  the  velocity  field,  upstream.  In  practical 
applications  to  part ia 1 ly-parahol ic  flows,  one  may  treat  the 
downstream  pressure  p^  as  a  known  quantity  with  =  0,  so  that  a 
marching,  instead  of  an  iterative  process,  can  he  employed  for 
the  calculation  of  the  transport  quantities.  The  pressure  field 
is  then  updated  by  repeating  the  marching  process  from  upstream 
to  downstream  for  several  sweeps  untiL  the  converged  pressure  and 
velocity  fields  are  obtained. 

In  the  present  study,  the  systems  of  algebraic  equations 
formed  by  the  assembly  of  the  pressure  and  the  pressure- 
correction  equations,  (IV-35)  and  ( IV— 39) ,  respec t i ve ly ,  are 
solved  by  the  tridiagonal-matrix  algorithm  with  six  line-by-line 
internal  iterations.  The  finite-analytic  coefficients  ap,  a j , 
etc.,  were  updated  in  each  upstream  to  downstream  global  sweep, 
but  remained  the  same  during  the  internal  iterations. 

IV. 4  Boundary  Conditions 

With  the  part ia 1 ly-parabol ic  assumptions,  the  transport 
equations  for  <P  -  U,V,W,k  and  e  become  parabolic  in  the  pre¬ 
dominant-flow  direction  C ,  but  remain  elliptic  in  the  cross- 
sectional  n4-  plane.  Due  to  the  latter  feature,  it  is  necessary 
to  specify  the  boundary  conditions  on  all  boundaries  of  each  04- 
plane.  On  the  other  hand,  the  parabolic  behavior  in  the  ;  - 
direction  enables  us  to  employ  a  marching  technique  without 
specifying  the  downstream  boundary  conditions  for  <t> .  For  the 
f ul ly-el 1 iptic  pressure  field,  one  must  specify,  either 
explicitly  or  implicitly,  the  boundary  conditions  on  all 
boundaries  including  the  downstream  exit  plane,  so  that  the 
upstream  transmission  of  pressure  information  is  properly 
represented.  In  other  words,  one  needs  to  specify  only  the 
pressure  datum  in  one  numerical  cell  if  all  velocity  components 
ire  proscribed  on  the  boundary  surfaces.  Alternatively,  when  the 
pressure*  is  prescribed  on  part  of  the  boundary  surfaces,  then 
inflow  or  outflow  must  be  permitted  there  to  avoid  overspec¬ 
ification  of  the  boundary  conditions. 


The  appropriate  boundary  conditions  for  ship-like  bodies  are 
as  f ol Lows . 

(1)  Initial  or  Upstream  Section  (4  =  1) 

The  distribution  of  the  velocity  components  (U,V,W)  and  the 
turbulence  quantities  (k  ,c)  are  prescribed  at  this  section  either 
from  detaiLed  thin  boundary- layer  calculations  or  from  simple 
flat-plate  correlations.  Boundary  conditions  for  pressure  and 
pressure-correction  are  not  required  in  the  present  staggered- 
grid  arrangement  since  the  pressure  is  implicitly  determined  by 
the  specification  of  the  velocity  components.  Also,  since  the 
velocity  components  are  known,  there  are  no  corrections  to  be 
made,  i  .  e  .  ( (.! ,  V  ,  W )  =  (  U  ,  V  ,  w  ) . ,  =  (  U  ,  V ,  W )  u  . 

(2)  K.xit  or  Downstream  Plane 

The  exit  plane  is  locate  :  in  the  wake  far  downstream  from 
the  stern,  and  the  zem  pressure-gradient  (p^  =  0)  condition  is 
specified  there.  Due  to  the  parabolicity  of  the  transport 

equations,  exit  conditio,  >.  for  <f>  =  U,V,W,k,  e  are  not  required. 

(3)  Body  Surface  (n=i> 

Per  laminar  flow,  the  numerical  solutions  are  carried  out 
upto  the  solid  surface  where  the  usual  no-slip  conditions, 

1 1  =  V  =  W - 0 ,  are  imposed. 

Kor  turbulent  flows,  a  modified  "wall-function"  approach 
(('hen,  Yoon  and  Yu  (1982))  is  employed  to  avoid  the  use  of  low- 
Reynold-number  turbulence  models  and  a  large  number  of  grid 
points  to  resolve  the  large  gradients  in  the  near-wall  region. 
The  boundary  conditions  are  determined  in  the  logarithmic  layer 
by  using  the  universal  1  a w-of - 1 he-wa 1 1  in  the  form 


U~  -  7  (ln  1 1-  11  -  ‘tYJi/2 - -1  *  21(1  +  4t/)1/2-  1]) 

t  t  (1  +  Ary  )  +  1 

+  B  +  3.7  A 


( IV-40 ) 
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where  U  is  the  friction  velocity  defined  by  u  =  /x  /pu  ' 

X  J  TWO 

y  =  Re  U  y  is  the  dimensionless  distance  measured  in  the 

T  3 

direction  normal  to  the  surface,  =  Vp/ReUT  is  the 

dimensionless  pressure  gradient,  A^  is  the  dimensionless  shear- 

stress  gradient  and  is  assumed  to  be  1/2  V  ,  q  is  the  magnitude 

P 

of  the  velocity,  <  =  0.42  is  the  von  Karraan  constant  and  B  = 

5.45.  In  the  present  study,  it  is  assumed  that  at  least  two 
points  ( n  =  2  and  3  in  Fig.  7)  are  located  in  the  logarithmic 
region.  This  avoids  the  need  for  a  separate  Couette-flow  type 
analysis  for  the  flow  between  the  wall  and  the  first  near-wall 
mesh  point  which  has  been  used  in  almost  all  previous 
applications  of  wall  functions. 

In  the  present  procedure,  a  value  for  is  assumed  and  the 
boundary  conditions  at  n  =  2  are  determined  from  the  logarithmic 
law  and  the  assumptions  of  local  equilibrium,  i.e., 


2  1  n  r4 

TT  =  7  {ln  [a 

T 


a  +  \y+2)1/2-  1  +  1/2 

- -  +  i/2 — i  +  +  ATy?>  -  D> 

X  (1  +  hxY2)  /  +  1 


+  B  +  3.7  A 


( IV-4 la ) 


(  IV-4 1 b ) 


( IV-4 lc ) 
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The  numerical  solution  then  provides  the  velocity  at  n  =  3  and 

IJ  is  updated  by  requiring  this  velocity  to  satisfy  the  log¬ 
arithmic  law,  i.e.,  by  Solving 
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+  B  +  3 . 7  A  ( TV-42) 

P 

by  a  root-finding  technique.  Thus,  an  iterative  procedure  is 
used  to  satisfy  the  wall  boundary  conditions  in  the  case  of  tur¬ 
bulent  flow,  and  to  determine  n  .  Typically,  five  iterations  are 
required  to  obtain  satisfactory  c  -nveryence.  In  the  present 
study,  the  pressure-gradient  correction  is  used  only  in  regions 
of  adverse  pressure  gradients. 

By  anchoring  the  solution  at  two  near-wall  mesh  points  on 

the  logarithmic  law,  the  present  wal 1 -f unct i on  approach  removes 
much  of  the  sensitivity  or  the  numerical  solution  to  the  location 
of  the  first  mesh  po  ■  ,t  which  has  been  observed  in  previous 
treatments.  The  pr  icedure  is  quite  straightforward  for  two- 

dimensional  and  a* isymmetr ic  flows.  For  three-dimensional  flows, 
however,  an  additional  assumption  concerning  the  direction  of  the 
velocity  vector  is  required  to  determine  the  individual  compo¬ 
nents  (U,V,W)  because  the  logarithmic  law  gives  only  the  velocity 

magnitude.  The  velocity  components  (U,V,W)  occurring  in  the 

£  n  C 

governing  equations  are  related  to  the  components  (q  ,  q  ,  q  ) 
along  the  body-fitted  coordinates  (Fig.  8)  by 

u  =  (U,V,W)  -  q  t,+  qn  Tr)+  q4  (IV-43) 


where,  from  equation  (A-10a), 
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and  the  inverse  relations  are 
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( IV- 4  4  ) 


(  TV-45) 


(  IV- 4 6  ) 


Now,  we  impose  the  contraint  that  q2  is  parallel  to  the  wall  and 

n 

therefore  require  q2  =  0  with  the  geometric  coefficients  deter¬ 

mined  at  n  =  1.  It  is  further  assumed  that  there  is  no  rotation 
of  the  velocity  vector  between  n  =  2  and  n  =  3  in  the  plane 

parallel  to  the  surface,  i.e. 
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(  IV- 4 7  ) 
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With  an  assumed  H  ,  and  therefore  q2  ,  equation  (IV-47)  gives 
C  i  Tfl 

and  q^  ,  and  with  q^  =  0,  equation  (IV-45)  gives  the  compo¬ 
nents  ( U  2  /  V  2  '  W2  '  required  as  the  boundary  conditions  in  the 
numerical  solution.  In  turn,  this  solution  gives  (U-^,V^,W^)  and 
hence  q^,  and  a  new  value  of  U  from  equation  (IV-42),  new  values 
of  q^  and  q^  from  equation  (IV-46),  and  finally  q^  and  q^  from 
equation  (IV-47)  which  provide  the  updated  boundary  conditions. 
This  procedure  is  iterated  to  cnnve.gence. 


Tt  should  also  be  noted  t  hat  in  evaluating  the  distance  y 

normal  to  the  wall,  require,,  in  the  law  of  the  wall,  it  is  neces- 

n 

sary  to  determine  the  direction  cosine  between  x^  and  n  .  This 
is  given  by 


COS  Cl 


2  2  2 

b. x„  +  b„r  +  b,r9„ 
In  2  n  3  n 
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2  2 


(  IV-48  ) 


with  all  the  geometric  coefficients  evaluated  at  n  =  1 . 

(4)  Outer  Roundary 

In  order  to  avoid  separate  viscous  and  inviscid  calcula¬ 
tions,  and  matching  between  them,  the  outer  boundary  is  placed  at 
a  large  distance  from  the  body  where  the  flow  is  assumed  to  be 
uniform.  The  boundary  conditions  then  become 

3k  3  e 

U=l,  W  =  0,  17=T^=0'  p=0  ( IV-49 ) 

The  radial  outflow  V  at  the  outer  boundary  is  not  specified  but 
is  determine'!  by  the  solution. 


*-  - 


-  v  ^ 


(5)  Symmetry  Planes 


On  the  ship  centerplane  and  on  the  undisturbed  free  surface 
or  waterplane,  symmetry  conditions  for  transport  quantities, 

i  .  e .  , 


W  =  0, 


3U  3y  _  3k  _  3e 
3?  "  3?  "  3C  "  3? 


are  enforced. 


(  IV- 50  ) 


(6)  Wake  Centerline 

For  axisymmetric  bodies,  the  boundary  conditions  at  the  wake 
centerline  r  =  0  (n  =  1)  are  specified  by 


V  =  W  =  0, 


3 U  3k  3  e 
3n  ~  3  n  ~  3 n 


(  IV- 51  ) 


For  three-dimensional  bodies,  the  above  zero-gradient 
conditions  are  no  longer  adequate  since  the  variations  of  U,  k 
and  e  in  the  circumferential  direction,  in  general,  result  in 
multiple  values  at  r  =  0.  In  order  to  avoid  this  problem,  we 
exclude  the  wake  centerline  in  the  solutions  and  apply  the  zero- 
gradient  boundary  conditions  on  the  surface  n  =  2. 


(7)  Initial  Conditions  (t  =  0) 

For  the  steady-flow  calculations  considered  here,  the 
initial  conditions  for  the  transport  quantities  <j>  =  U,V,W,k, 
are  taken  directly  from  the  values  known  at  the  immediately  up¬ 
stream  station  during  the  first  sweep  and  the  pressure  is  assumed 
to  be  zero  throughout  the  flow  field.  Although  it  is  possible  to 
assume  more  realistic  initial  conditions  and  thus  accelerate  the 
convergence  of  the  solution,  these  rather  crude  initial  condi¬ 
tions  have  been  used  not  only  to  simplify  the  use  of  the  method 
but  also  to  demonstrate  its  versatility. 


IV. 5  The  Overall  Solution  Algorithm 

The  solution  algorithm  presented  below  differs  from  that 
outlined  in  Chen  and  Patel  (  1984)  in  the  way  the  pressure  is 
updated.  Tn  the  previous  procedure,  the  pressure  field  was  ob¬ 
tained  station  by  station  by  solving  the  momentum  and  pressure 
equations  simultaneously  during  the  marching  process.  While  this 
gave  quite  sa t i s f actot y  results,  it  allowed  only  one  upstream-to- 
downstream  sweep  for  the  pressure  field  in  each  global  sweep,  and 
therefore  the  pressure  influence  propagated  only  one  station 
upstream  per  global  sweep.  The  number  of  sweeps  required  to 
achieve  convergence  with  this  procedure  was  typically  of  the 
order  of  the  number  of  streamwise  stations  employed  in  the  cal¬ 
culations. 

Tn  order  to  improve  the;  fate  of  convergence,  a  new  global 
pressure- i terat ion  technique  is  employed  here  to  bring  the  down¬ 
stream  pressure  influence  immediately  to  the  whole  solution  do¬ 
main.  The  new  solution  algorithm  requires  only  a  simple  rear¬ 
rangement  of  the  solution  procedures  of  Chen  and  Patel  (1984)  and 
a  slight  increase  of  computer  storage  for  pressure  coeffi¬ 
cients.  In  the  new  global  pressure- iterat ion  scheme,  the  momen¬ 
tum  equations  are  solved  using  the  old  pressure  field  during  each 
global  sweep.  After  the  sweep  for  the  transport  quantities  is 
completed,  the  pressure  field  is  solved  elliptically  with  several 
global  iterations  just  as  in  f ul ly-el 1  ipt ic  problems.  These 
iterations  can  be  performed  from  downstream  to  upstream  so  that 
the  downstream  pressure  influence  propagates  all  the  way  to  the 
inlet  station  even  in  the  first  global  sweep.  In  this  fashion, 
the  number  of  global  sweeps  required  to  achieve  the  same  level  of 
convergence  as  before  is  reduced  by  a  factor  of  two  or  more, 
while  the  computation  time  for  each  global  sweep  remains  about 
the  same.  Hence,  the  overall  computation  time  is  greatly  re¬ 
duced. 

For  transient  problems,  where  the  initial  and  boundary  con¬ 
ditions  are  properly  specified  as  functions  of  time,  the  overall 
numerical  solution  procedure  may  be  summarized  as  follows: 


Construct  the  body-fitted  coordinate  system  for  the  given 
body  shape  and  solution  domain,  and  calculate  the  geometric 
coefficients  b|  ,  g1^,  J,  etc.  from  equations  (IV-27)  to  (IV- 
31  )  . 

Specify  the  initial  conditions  for  the  velocity  and  turbu¬ 
lence  fields.  Set  p  =  0  everywhere  initially. 

Specify  the  velocity  and  turbulence  profiles  at  the  first 
station  =  1  (these  may  be  time  dependent). 

Calculate  the  finite-analytic  coefficients  for  pressure, 
pressure-correction  and  momentum  equations  at  the  downstream 
station.  Store  only  the  finite-analytic  coefficients  a^  ,  an 
and  a  for  the  pressure  equation. 

Solve  the  momentum  equations  (IV-29)  based  on  the  updated 
pressure  field  to  obtain  the  starred  velocity  field  (U*,  V* , 
V/* )  .  This  system  of  algebraic  equations  is  solved  by  a  tri¬ 
diagonal  matrix  algorithm. 

Calculate  the  mass  source  D* ,  and  solve  the  pressure  correc¬ 
tion  equation  (IV-39)  by  tridiagonal  matrix  algorithm. 

Correct  the  velocity  field  using  the  velocity-correction 
formulae  (IV-38),  but  do  not  correct  the  pressure  field. 
Calculate  the  pseudoveloc ities  (u,  V,  W)  in  terms  of  the 
velocity  field  from  equation  (IV-33).  Store  only  D  for  later 
use . 

Solve  equations  (IV-29)  for  turbulence  quantities  ( <t>  =  k,e) 
by  tridiagonal-matrix  algorithm. 

March  to  the  next  downstream  station  and  repeat  steps  (3)  to 
(9)  . 

After  reaching  the  last  downstream  station,  solve  the  pres¬ 
sure  equation  (IV-35)  by  tridiagnonal  matrix  algorithm. 
Several  iterations  from  downstream  to  upstream  are  employed 
to  update  the  three-dimensional  elliptic  pressure  field. 
Repeat  steps  (3)  to  (11)  for  several  sweeps  until  both  the 
pressure  and  velocity  fields  have  converged  within  a  speci¬ 
fied  tolerance. 


(13)  Return  to  step  (2)  for  the  next  time  step. 

(14)  Stop  if  the  steady-state  solution  is  achieved,  or  if  time 
exceeds  the  maximum  time  period  assigned.  For  steady-flow 
calculations,  one  may  relax  the  convergence  criterion  in  step 
(12)  and  use  a  larger  time  increment  for  the  intermediate 
solutions. 

In  the  present  study,  only  one  sweep  was  used  in  step  (12) 
for  each  time  step.  Also,  as  noted  earlier,  instead  of 
specifying  the  initial  conditions  for  velocity  and  turbulence 
profiles  everywhere  at  t  =  0,  one  may  specify  only  the  profiles 
at  the  first  station  5=1.  The  downstream  profiles  are  then 
taken  from  the  immediate  upstream  station  (i.e.,  =  4>  at  t  = 

0)  during  the  first  global  sweep. 
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V.  SOME  EXAMPLE  SOLUTIONS 

The  numerical  method  described  in  Section  IV  has  been 
employed  to  calculate  a  variety  of  two-dimensional,  axisymmetric 
and  three-dimensional  flows.  Here,  we  shall  present  only  a  few 
representative  examples  which  illustrate  the  capability  of  the 
method  and  indicate  areas  which  require  further  study. 

V. 1  Two-Dimensional  Flow  on  a  Flat  Plate 

This  is  by  far  the  simplest  example  of  flows  which  involve 
the  evolution  of  a  wake  from  a  boundary  layer  because  the  vis- 
cous-inviscid  interaction  at  the  trailing  edge  is  weak,  the  flow 
is  two-dimensional,  and  Cartesian  coordinates  can  be  used  in  the 
physical  plane.  Nevertheless,  it  provides  a  useful  test  of  the 
grid-generation  technique,  the  numerical  method,  and  the  turbu¬ 
lence  model. 

(i)  Laminar  Flow 

The  laminar  flow  over  the  trailing  edge  of  a  flat  plate  has 
been  the  subject  of  numerous  previous  investigations  using  a 
variety  of  methods.  Among  these  are  numerical  solutions  of  the 
full  Nav ier-Stokes  equations  (Janssen,  1958;  Dennis  and  Dunwoody, 
1966)  for  relatively  low  Reynolds  numbers,  numerical  solutions  of 
the  triple-deck  equations  of  Stewartson  (1969)  and  Messiter 
(1970)  which  apply  at  very  large  Reynolds  numbers  (Jobe  and  Burg- 
yraf,  1974;  Melnik  and  Chow,  1975;  Veldman,  1975),  numerical 
solutions  of  the  part ia 1 ly-parabol ic  equations  (Rubin  and  Reddy, 
1983;  Saint-Victor  and  Cousteix,  1984),  and  solutions  employing 
viscous-inviscid  interaction  techniques  in  which  the  boundary- 
layer  equations  are  coupled  with  the  external  inviscid  stream 
(Veldman,  1979,  1981;  Davis  and  Werle,  1981).  A  thorough  review 
of  the  similarities  and  differences  among  the  various  approaches 
and  their  results  is  beyond  the  scope  of  the  present  paper. 
However,  the  results  of  the  present  method  will  be  compared  with 
some  representative  solutions  obtained  earlier. 


We  select  Cartesian  coordinates  (x,y)  in  the  physical  plane, 
with  the  origin  at  the  leading  edge.  With  distances  normalized 
by  the  pi  a  length  I.,  x  =  0  and  I  correspond  to  the  leading  and 
trailing  edges,  respectively,  and  y  is  normal  to  the  plate.  A 
nonuniform  rectangular  gria  was  generated  using  the  technique 
described  in  Section  IV. 1  with  x  =  x(£)  and  y  =  y ( n )  .  With  this 
rather  simple  coordinate  transformation,  equations  ( IV— 6 )  simpli¬ 
fy  to 


1  1 


f  x  .  =  0 
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(V-l  ) 


If  the  required  grid  distribution  in  the  physical  space  is  chosen 
a  priori,  then  the  control  functions  can  be  determined  from 
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( V-  2 ) 


Alternatively,  one  may  prescribe  the  control  functions  and  solve 
equations  (V-l)  for  the  coordinates  (x,y).  For  the  calculations 
presented  here,  the  latter  approach  is  adopted  to  ensure  a  smooth 
variation  of  the  coordinates  and  the  control  functions.  The 
particular  functions  used  are 
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where 


5  -  25j  +  1 
z2  =  5  -  25  +  1 


5^  and  correspond  to  the  leading  and  trailing  edges,  respec¬ 
tively,  and  A-j^  ,  A2  and  A^  are  positive  constants  which  can  be 
selected  to  yield  the  desired  grid  concentration  around  x  =  0,  1 
and  y  =  0. 


If  we  further  assume  that  the  grid-control  functions  remain 
constant  within  each  numerical  cell  shown  in  Fig.  9,  equations 
( V— 1 )  can  be  solved  analytically  to  obtain 


a  t  -a 
e  XU  +  e  *D 
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( V-4  ) 


where 
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b 
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Equations  (V-4)  constitute  a  set  of  simultaneous  algebraic  equa¬ 
tions  which  can  be  solved  by  the  tridiagonal-matrix  algorithm. 
The  first  derivatives  of  the  coordinates,  which  appear  in  the 
geometric  coefficients,  g1  -1  ,  b|  ,  etc.,  can  then  be  obtained 
readily  from 
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The  numerical  grid  shown  in  Fig.  10  was  generated  with 

Aj  =  0.3,  A2  =  0*2,  A^O.4,  and  5  =  1  at  x  =  -1.385,  5=19  ( =  5^)  at 

x  =  0,  5  =  49  (  =  £,  )  af;  X  =  1  >  5  =  64  at  x  =  3  .  339  ,  n  =  2  at  y  =  0  and  n=16  at 

y=12.  Thus,  a  64  x  15  mesh  covers  the  physical  region  that 


t 


< 
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extends  from  a  distance  1.385  L  upstream  of  the  leading  edge  to 
2.339  L  downstream  of  the  trailing  edge,  and  12  L  normal  to  the 
plate,  with  the  grid  concentrated  in  the  neighborhood  of  the 
leading  and  trailing  edges  and  the  plate.  The  domain  can  be 
readily  extended  in  all  directions  by  the  addition  of  only  a  few 
extra  nodes. 

The  present  method  has  been  employed  with  different  grid 
distributions  and  mesh  sizes,  and  different  locations  of  the 
upstream,  downstream  and  outer  boundaries  to  study  their  influ¬ 
ence  on  the  solution.  These  will  be  discussed  in  a  separate 
paper.  For  the  present  purposes  it  suffices  to  present  and  com¬ 
pare  solutions  obtained  wirh  grid  distributions  and  solution 
domains  similar  to  those  employed  by  others.  Thus,  we  shall 
concentrate  on  the  flow  over  the  trailing  edge  and  employ  only  a 
portion  of  the  grid  described  above. 

As  in  some  previous  studies,  the  par t i a  1 ly-parabol ic  equa- 

5 

tions  have  been  solved  for  a  plate  Reynolds  number  UoL/v  =  10  in 
the  domain  0.4819  <  x  <  1  . 4 3 ,  0  <  y  <  0.2196  (i.e.  37  <  5  <  60, 

2  <  n  <  11),  which  c  •  esponds  roughly  to  the  recent 
calculations  of  Saint-Victor  and  Cousteix.  At  the  upstream  sta¬ 
tion,.  the  streamwise  component  of  velocity  t.)  is  specified  by  the 
Blasius  solution  and  the  streamwise  gradient  of  the  normal  com¬ 
ponent,  Vx ,  is  set  equal  to  zero.  All  other  boundary  conditions 
are  prescribed  as  discussed  in  Section  IV. 4. 

Figure  11  shows  the  pressure  distribution  on  the  plate  and 
along  the  wake  centerline  calculated  at  different  time  steps  (or 
global  iterations),  starting  with  an  initially  zero  pressure 
throughout.  It  is  seen  that  the  solution  converges  monoton ica 1 ly 
in  about  20  steps,  the  convergence  being  more  rapid  over  the 
plate  than  in  the  wake.  The  rapid  monotonic  convergence  of  the 
solution  is  believed  to  .be  due  to  the  two-step,  global  pressure- 
correct  ion  procedure  described  in  Section  IV. 3.  Sol ut ions  with 
some  previous  algorithms  have  led  to  an  oscillatory  behavior  of 
pressure,  requiring  several  hundred  iterations  for  acceptable 
converge nee . 


The  important  features  of  the  present  solution  are  shown  in 
Figs.  12,  13  and  14,  along  with  the  interacting  boundary- layer 
calculations  of  Veldman  (1979),  the  triple-deck  results  of  Melnik 
and  Chow  (1975)  and  the  pact ial ly-parabol ic  calculations  of 
Saint-Victor  and  Cousteix  (1984).  Although  all  four  solutions 
are  qualitatively  similar,  a  closer  examination  reveals  several 
interesting  differences. 

First,  we  consider  the  flow  on  the  plate  (x  <  1).  It  is 
seen  that  the  triple-deck  calculations  predict  a  lower  pressure 
than  that  calculated  by  the  other  methods  over  the  middle  of  the 
plate.  This  is  due  to  the  fact  that  the  triple-deck  solutions 
are  obtained  with  the  Blasius  solution  prescribed  at  an  asympto¬ 
tically  large  distance  upstream  of  the  trailing  edge,  whereas  the 
other  three  methods  specify  the  initial  conditions  at  x  ~  0.5. 

5 

In  other  words,  for  the  Reynolds  number  of  10  considered  here, 
the  triple-deck  results  suggest  that  the  influence  of  the 
trailing  edge  penetrates  further  upstream.  In  fact,  other 
calculations  performed  with  the  present  method  with  initial  con¬ 
ditions  closer  to  the  leading  edge  (typically  x  =  0.18)  lead  to 
results  in  much  better  agreement  with  the  triple-deck  solutions, 
and  considerably  reduce  the  difference  between  the  present  solu¬ 
tions  and  those  of  Veldman,  and  Saint-Victor  and  Cousteix,  at  the 
trailing  edge.  The  latter  authors  also  commented  on  the  in¬ 
fluence  of  the  location  of  the  outer  boundary  of  the  solution 
domain  on  the  upstream  pressure  distribution.  The  present  method 
also  indicated  such  a  sensitivity.  The  wall  shear-stress  coef¬ 
ficient  (C  =  r  Re  t  /'otl  )  shown  in  Fjq.  II  Indicates  somewhat 
TWO 

better  agreement  in  the  results  >.f  all  methods  except  that  the 
present  solution  predicts  higher  values  close  to  the  trailing 
edge,  as  might  he  expected  from  the  larger  favorable  pressure 
gradient  calculated  by  this  method. 

Turning  next:  to  the  flow  in  the  no  a  r  wake,  it  is  seen  that 
the  pressure  distribution  obtained  with  the  present:  method  agrees 
rather  wo*  1  1  with  that,  of  sa  i  nt -v  i  ct.or  and  Cousteix.  but  both 
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predict  higher  pressures  than  the  triple-deck  and  interacting 
boundary-layer  theories.  On  the  other  hand,  the  velocity  varia¬ 
tion  along  the  wake  centerline  (Fig.  14)  shows  quite  different 
trends.  Here,  the  solutions  of  Saint-Victor  and  Cousteix  agree 
with  those  of  Veldman,  the  present  calculation  predicts  higher 
velocities,  while  the  triple-deck  solution  diverges  from  all 
others.  The  latter  feature  is  obviously  due  to  the  fact  that  the 
triple-deck  solution  is  matched  downstream  with  the  Gold¬ 
stein  near-wake  rather  than  the  x^/^  parabolic  asymptotic  wake. 
The  remaining  solutions  reproduce  the  parabolic  characteristic 
due  to  the  downstream  boundary  condition  px  =  0  imposed.  How¬ 
ever,  the  reasons  for  the  disagreement  in  wake  velocities  inspite 
of  the  agreement  in  the  pressure  distributions  of  the  present 
method  and  that  of  Saint-Victor  and  Cousteix,  and  for  the  disa¬ 
greement  in  the  pressure  distribution  inspite  of  the  agreement  in 
the  wake  velocity  of  the  methods  of  Veldman  and  Saint-Victor  and 
Cousteix,  are  not  clear.  Calculations  performed  with  the  present 
method  by  moving  the  downstream  and  outer  boundaries  further  away 
from  the  plate  (x  =  3.34  and  y  =  12.0,  respectively)  bring  the 
predicted  wake  pressure  distribution  into  better  agreement  with 
the  triple-deck  and  Veldman  solutions  but  the  wake  velocity 
distribution  remains  unaffected.  This  leads  us  to  conclude  that 
the  solutions  of  the  part ially-parabol ic  equations,  both  present 
and  past,  are  sensitive  to  the  location  of  the  solution 
boundaries  and  agreement  with  the  triple-deck  solution  can  be 
secured  only  by  employing  a  large  enough  solution  domain.  These 
aspects  will  be  discussed  in  a  separate  report. 

(ii)  Turbulent  flow 

The  calculations  for  the  turbulent  flow  over  the  trailing 
edge  of  a  flat  plate  presented  below  correspond  to  the  experi¬ 
ments  of  Ramaprian,  Patel  and  Sastry  (1981).  If  X  and  Y  are 
coordinates  measured  along  and  normal  to  the  plate  with  X  =  0 
being  the  trailing  edge,  L  (=  1.829  m)  is  the  length  of  the 


plate,  and  U0  (=  22  m/s)  is  the  freestream  velocity,  the  calcula¬ 
tions  correspond  to  a  Reynolds  number  UqL/v  =  1.6  x  10^. 

In  order  to  obtain  the  desired  grid  concentration  in  the 
neighborhood  of  the  trailing  edge  and  the  plate,  the  numerical 
grid-generation  technique  was  again  employed  but  with  and  C2 
corresponding  to  the  initial  station  in  the  boundary  layer  at  X/L 
=  -  0.6,  and  the  trailing  edge  X  =  0,  respectively.  The  grid 
used  here  corresponds  to  A1=0.20,  A2=0.12  and  A3=0.25; 

with  £  =  «i1  =  7  at  X=-0 .6  L,  5  =  =31  at  X=0,  5  =63  at  X=8.57 

L,  n  =1  at  Y=0  and  n  =15  at  Y=L.  In  other  words,  a  57  x  15  mesh 
is  used  to  cover  a  physical  region  of  -  1097  mm  <  X  <  15680  mm 
(-0.6  <X/L  <8.57;  7  <  £  <63)  and  0  <Y  <  1829  mm  ( 0  <  Y/L  <  1.0; 

1  <  n  <  15).  Standard  flat-plate  boundary-layer  profiles  were 

specified  at  the  upstream  station. 

Since  we  are  interested  only  in  the  steady-state  solution,  a 
very  large  time  increment,  t  =  100,  was  used  with  only  one  global 
sweep  per  time  step.  The  convergence  histories  of  several  quan¬ 
tities  of  particular  interest  are  shown  in  Figs.  15  through  18  at 
nine  different  time  steps  corresponding  to  sweep  numbers  1,  10, 

20,  30,  40,  50,  60,  80,  and  100.  It  can  be  seen  that  all  quan¬ 
tities  converge  monotonically  in  about  30  sweeps  although  the 
calculations  were  continued  for  100  sweeps  or  time  steps. 

Figure  15  shows  the  pressure  distribution  on  the  plate  and 
along  the  wake  centerline.  It  is  seen  that  a  sharp  reduction 
followed  by  a  rapid  increase  of  pressure  is  predicted  near  the 
trailing  edge  as  a  result  of  the  sudden  change  fr<->m  the  no-slip 
boundary  condition  on  the  plate  to  the  zero-shear-stress  cond¬ 
ition  along  the  wake  centerline.  Comparison  of  Fig.  15  with  the 
corresponding  case  for  laminar  flow,  Fig.  12,  indicates  that  the 
viscous-inviscid  interaction  in  turbulent  flow  is  stronger,  as 
measured  by  the  pressure  changes,  but  confined  to  a  smaller 

region  around  the  trailing  edge.  The  wall  shear  stress  coef- 

.  .  1  2 

ficxent  (Cf  =  i  /  PUQ) ,  shown  in  Fig.  16,  decreases  continously 
toward  the  trailing  edge  during  the  first  parabolic  sweep  under 
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zero  initial  pressure  field.  However,  subsequent  iterations  lead 

to  an  increase  near  the  trailing  edge  as  in  the  laminar  case. 

The  higher  shear  stress  indicates  an  acceleration  of  the  fluid 

inside  the  boundary  layer  ahead  of  the  trailing  edge  and  is 

associated  with  the  favorable  pressure  gradient  in  this  region. 

It  is  interesting  to  note  that  both  the  turbulent  kinetic  energy 

(k)  and  the  rate  of  turbulent  energy  dissipation  (e),  and 

therefore  the  centerline  eddy  viscosity  shown  in  Fig.  18, 

decrease  considerably  during  the  time-marching  process.  This  is 

obviously  related  to  the  decrease  of  velocity  gradients  in  the 

normal  direction  as  the  solution  converges.  When  convergence  is 

2 

achieved,  the  eddy  viscosity  { =  c^k  /e)  along  the  wake 

centerline  increases  rapidly  with  distance  from  the  trailing  edge 

-4 

and  reaches  an  essentially  constant  value  of  v  /u  l  =  1.74  x  10 
beyond  X  >  6000  mm  (i.e.,  X/L  >  i .  1 )  .  The  significance  of  this 
result  is  discussed  later  on. 

The  principal  results  of  the  calculations  in  the  near  wake 
region  are  compared  with  experimental  data,  which  extended  only 
upto  X  =  610  mm,  in  Figs.  19,  20  and  21.  The  level  of  agreement 
between  the  calculation  .  and  experiments  with  respect  to  the 
increase  of  centerline  elocity  (Fig.  19),  mean-velocity  profiles 
(Fig.  20),  and  turbulent  kinetic-energy  distributions  (Fig,  21) 
is  comparable  with  that  achieved  by  Patel  and  Scheuerer  (1982), 
who  solved  the  boundary- layer  equations  using  the  same  turbulence 
model  but  a  different  numerical  method.  Computationally,  the  two 
methods  differ  in  the  number  of  grid  points  required  to  resolve 
the  large  gradients  in  the  near  wake.  The  solutions  of  Patel  and 
Scheuerer  were  obtained  with  60  points  across  the  wake  and  some 
800  streamwise  steps  in  the  region  0  <  X  <  610  mm  whereas  the 
present  solutions,  as  noted  above,  employed  57  x  15  grid  nodes 
over  a  much  larger  domain.  In  fact,  as  shown  in  Fig.  21,  only  8 
points  in  the  y-direction  are  needed  to  resolve  the  boundary 
layer  and  the  near  half-wake,  whose  thickness  is  of  the  order  of 
50  mm.  The  remaining  7  points  cover  the  inviscid  domain  with  the 
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outer  boundary  located  at  Y  =  L  =  1829  mm.  This  difference  be¬ 
tween  the  nodes  required  in  the  two  numerical  schemes  is  at¬ 
tributed  largely  to  the  use  of  the  finite-analytic  scheme  in  the 
present  method. 

Physically,  the  two  methods  are  also  different.  The  use  of 
the  boundary-layer  equations  neglects  the  viscous-inviscid  inter¬ 
action  at  the  trailing  edge  while  the  partially-parabol ic  equa¬ 
tions  take  it  into  account.  In  the  application  of  boundary-layer 
equations  to  the  flow  on  a  flat  plate,  it  is  assumed  that  the 
pressure  is  uniform  (equal  to  zero,  say)  everywhere  and, 
therefore,  the  velocity  outside  the  boundary  layer  and  wake  is 
constant  and  equal  to  the  freestream  velocity  UQ.  In  the  par¬ 
tially-parabol  ic  solutions  the  pressure  is  set  equal  to  zero  only 
at  very  large  distances  from  the  plate  and,  associated  with  the 
predicted  pressure  distribution  on  the  plate  and  along  the  wake 
centerline,  shown  in  Fig.  15,  is  a  weak  but  nonzero  pressure 
field  everywhere.  In  the  inviscid  nonturbulent  flow  outside  the 
boundary  layer  and  wake,  this  results  in  velocities  greater  than 
U0,  i.e.,  an  "overshoot"  in  the  velocity  profiles.  Although  the 
overshoots  are  quite  small  in  magnitude,  the  maximum  being  of  the 
order  of  0.007  UQ  in  the  present  case  (these  are  barely  visible 
in  traditional  plots  like  Fig.  20),  and  decrease  to  zero  with 
distance  from  the  trailing  edge  in  all  directions,  they  are  of 
considerable  significance  in  the  interpretation  of  the  results  in 
the  far  wake  and  in  comparisons  with  the  so-called  asymptotic 
laws. 

The  overshoot  in  the  velocity  profiles  mentioned  above  im¬ 
mediately  lead  to  problems  in  the  definition  of  parameters  used 
in  the  traditional  analysis  of  asymptotic  wakes.  These  are  il¬ 
lustrated  in  Fig.  22.  The  assumption  of  constant  velocity  out¬ 
side  the  wake,  adopted  in  boundary-layer  theory  and  assumed  in 
experiments,  yields  the  usual  definitions  of  the  centerline  velo¬ 
city  defect  WQ  (=  UQ  —  Uc)  and  wake  half-width  b  (=  Y  where  W/WQ 
=  1/2).  These  are  shown  in  Fig.  22(a).  However,  if  the 


pressure,  and  therefore  the  velocity,  varies  in  the  normal 
direction  outside  the  wake,  then  the  physically  most  meaningful 
definitions  are  those  shown  in  Fig.  22(b).  In  this  case  the 
velocity  defect  is  the  difference  between  the  actual  velocity  and 
the  inviscid-f low  velocity  extrapolated  into  the  wake.  Note  that 
it  is  this  velocity  defect  which  is  related  to  the  drag 
coefficient  of  the  plate.  The  determination  of  the  so-defined 
deficit  from  a  numerical  solution  is  not  a  straightforward  matter 
because,  as  we  have  already  noted,  the  overshoot  is  small  in 
absolute  terms  but  not  small  in  comparison  with  the  velocity 
defect  which  decreases  with  downstream  distance  and  is  a 
difference  between  two  large  quantities.  In  view  of  these 
difficulties,  we  shall  adopt  the  conventional  definitions  even  in 
the  presence  of  an  overshoot  and  use  the  velocity  defect  WQ  (=  UQ 
-  Uc )  ,  and  the  corresponding  half-width  b  ,  as  shown  in  Fig. 
22(c).  These  have  the  advantages  that  they  can  be  determined 
precisely,  can  be  used  to  assess  certain  asymptotic  features  of 
the  solutions,  and  also  enable  a  comparison  to  be  made  between 
the  calculations  and  experimental  data  processed  with  the  same 
def initions. 

Figures  23  and  24  show  the  calculated  velocity  and  Reynolds 

shear-stress  profiles  at  several  stations  far  downstream  of  the 

trailing  edge  (X  >  2.3  L)  in  the  usual  defect  forms,  i.e., 

★  * 

using  WQ  and  b  as  the  character istic  velocity  and  length  scales, 
respectively.  The  velocity  overshoot  in  the  outer  flow  is 
clearly  evident  in  Fig.  23  and  amounts  to  about  12  percent  of  the 
centerline  defect.  Consistant  with  the  eddy-viscosity  assumption 
and  the  velocity  overshoot,  there  is  also  a  region  of  very  small 
negative  stress,  but  the  stress  reduces  to  practically  zero  at 
Y/b*  ~  1.2.  A  particularly  noteworthy  feature  is  that  the  velo¬ 
city  profiles  become  self-similar  somewhat  earlier  than  the  Rey- 
nold-stress  profiles,  with  similarity  in  both  beyond  about  6000 
mm  (i.e.  X  >  3.3  L,  approximately).  This  coincides  roughly  with 
the  region  beyond  which  the  wake-centerline  eddy  viscosity  became 
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constant  (see  Fig.  18).  Such  a  similarity  is  to  be  expected  on 
theoretical  grounds  because  the  partially-parabol ic  equations 
also  reduce  to  the  asymptotic,  small-defect  wake  equations  which 
are  usually  obtained  from  the  boundary-layer  equations. 
Ramparian,  Patel  and  Sastry  (1982)  reviewed  the  available  flat- 
plate  wake  data  and  concluded  that  the  asymptotic  behavior  is 
realized  with  respect  to  the  mean  velocity  and  turbulence  at  a 
distance  of  the  order  of  3500  from  the  trailing  edge,  6  being  the 
wake  momentum  thickness.  The  value  of  3.3  L  quoted  above  for  the 
present  calculations  corresponds  to  about  7709  and,  therefore, 
the  present  solutions  indicate  a  somewhat  slower  approach  to  the 
asymptotic  state. 

Perhaps  the  most  surprising  result  to  emerge  from  the  pre¬ 
sent  calculations  is  the  asymptotic  value  of  the  eddy  viscosity 

in  the  far  wake.  The  calculated  value  of  v  /u  L  =  1.74  x  10-4, 

t  o 

quoted  earlier  (see  Fig.  18),  corresponds  to  Vt/Uo0  =  0.04, 

which  is  somewhat  higher  than  0.035  (±  0.001)  deduced  by  Patel 

and  Scheuerer  (1982)  from  the  data  of  Pot  (1979)  in  the  range 

400  <  X/6  <  1000.  This  difference  is  also  observed  in  the  maxi- 

_  *  o 

mum  Reynolds  stress:  the  calculated  maximum  of  -  uv  /W^  being 
0.054  while  the  measurements  of  Pot  indicate  0.050.  These  re¬ 
sults  are  surprising  because  the  extensive  calculations  of  Patel 
and  Scheuerer,  who  used  boundary-layer  equations  but  the  same 
turbulence  model,  led  to  much  smaller  values  of  the  eddy- 
viscosity  (0.024)  and  maximum  shear  stress  (0.034).  This  sub¬ 
stantial  underpred ict ion  was  attributed  by  Patel  and  Scheuerer  to 
a  deficiency  of  the  k-e  turbulence  model.  The  present  calcula¬ 
tions,  however,  suggest  that  the  disagreement  between  the  experi¬ 
ments  and  calculations  based  on  the  boundary- layer  equations  may 
not  be  entirely  due  to  the  turbulence  model.  Rather,  the  varia¬ 
tion  of  pressure  and  the  normal-stress  terms  which  are  incorpor¬ 
ated  in  the  part ial ly-parabol ic  equations  but  neglected  in  the 
boundary-layer  equations  may  be  responsible  for  the  dif¬ 
ferences.  The  uncertainties  in  the  definitions  of  the  velocity 


and  length  scales  used  in  the  asymptotic  analysis  may  also  con¬ 
tribute  to  the  observed  differences.  Needless  to  say,  the  resol¬ 
ution  of  these  rather  perplexing  differences  is  quite  important 
in  order  to  assess  the  generality  of  the  turbulence  model.  Addi¬ 
tional  calculations  and  more  detailed  comparisons  are  obviously 
needed.  These  are  in  progress. 

In  conclusion,  we  note  that  the  part ial ly-parabol ic  method 
successfully  predicts  all  of  the  features  of  laminar  and  turbu¬ 
lent  flows  over  a  flat  plate.  Inspite  of  the  relative  simplicity 
of  this  two-dimensional  problem,  it  has  provided  an  opportunity 
to  test  several  basic  components  of  the  method  by  comparing  the 
present  results  with  those  obtained  by  others  with  quite  diffe¬ 
rent  methods. 

V.2  Thick  Axisymmetric  Boundary  Layer  and  Wake 

Calculations  using  the  present  method  have  been  made  for 
five  different  bodies  for  which  detailed  experimental  data  are 
available.  We  shall,  however,  present  the  results  for  only  two 
cases,  namely  the  models  with  Afterbody  5  and  Afterbody  1  tested 
by  Huang  et  al.  (1980,  1979).  The  former  is  of  particular  inter¬ 
est  since  the  stern  contains  a  point  of  inflextion  and  leads  to 
quite  dramatic  changes  in  the  surface  pressure  distribution.  The 
latter  is  included  here  becasue  it  is  the  axisymmetric  parent 
model  of  a  family  of  three-dimensional  bodies  of  elliptic  cross- 
sections  to  be  discussed  later. 

The  calculations  for  Afterbody  5  were  performed  with  60  sta¬ 
tions  in  the  domain  0.60  <_  X/L  13.20,  where  X  is  measured  along 
the  axis  of  the  body  from  the  nose  and  L  (=  9.55  ft)  is  the  body 
length,  and  19  points  between  the  body  and  the  external  boundary 
at  R/L  =  0.8137,  R  being  the  radial  distance  from  the  body 
axis.  A  partial  view  of  the  body-fitted  coordinates  is  shown  in 
Fig.  4.  The  use  of  coordinate-stretching  functions  in  the  longi¬ 
tudinal  and  radial  directions  ensures  that  the  grid  points  are 
closely  spaced  inside  the  viscous  region  and  near  the  stern. 


These  calculations  were  started  by  specifying  the  grid  geo¬ 
metry,  and  the  velocity  and  turbulence-parameter  profiles  in  the 
boundary  layer  at  the  upstream  station,  and  assuming  that  a  con¬ 
stant  ambient  pressure,  i.e.  p  =  0,  prevails  throughout  the  solu¬ 
tion  domain.  For  simplicity,  flat-plate  boundary-layer  profiles 
were  used  at  the  upstream  station  X  =  0.6  L,  where  the  boundary 
layer  is  thin.  The  numerical  solutions  were  then  obtained  using 
the  time-marching  technique  with  a  time  increment  of  x  =  1.0. 

Since  the  influence  on  the  transport  quantities  from  down¬ 
stream  comes  solely  from  the  elliptic  pressure  field,  a  converged 
solution  is  obtained  only  when  the  pressure  field  converges.  It 
is  therefore  very  important  to  examine  first  the  convergence 
history  of  the  pressure  field.  Fig.  25  shows  the  pressure  dis¬ 
tribution,  Cp  =  2p,  on  the  body  surface  and  along  the  wake 
centerline  calculated  at  different  time  steps  or  global  sweeps. 
It  is  seen  that  the  solution  converges  monotonically  and  a 
"steady-state"  or  converged  solution  is  obtained  after  about  30 
time  steps  or  global  sweeps.  The  converged  solution  is  in  fairly 
good  agreement  with  the  data. 

Fig.  26  shows  the  variation  of  pressure  across  the  boundary 
layer  and  in  the  external  flow,  as  a  function  of  the  distance 
from  surface,  (R  -  RQ) ,  where  RQ  is  the  local  radius  of  the  body, 
all  the  way  to  the  external  boundary  of  the  solution  domain, 
where  Cp  =  0  is  specified.  Here  again  the  agreement  with  the 
available  experimental  data  is  quite  satisfactory,  considering 
the  difficulties  of  accurately  measuring  pressure  in  such  an 
environment.  Note  that  uniform  ambient  pressure  is  recovered 
beyond  a  distance  of  about  0.35  L,  from  the  body  surface,  and 
consequently  we  conclude  that  the  region  of  viscous-inviscid 
interaction  in  this  case  extends  to  about  0.35L  or  4  maximum  body 
diameters.  Since  the  domain  of  the  present  calculations  extends 
well  beyond  this,  the  viscous-inviscid  interaction  is  completely 
captured  without  separate  viscous-flow  and  potential-flow 
calculations,  and  iterative  matching  between  them. 
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Detailed  comparisons  between  the  calculated  and  experimental 
profiles  of  the  axial  (U)  and  radial  (V)  components  of  mean  velo¬ 
city,  normalized  by  UQ,  shown  in  Fig.  27  indicate  that  the  boun¬ 
dary-layer  thickness  is  correctly  predicted.  We  also  note  that 
both  components  of  velocity  continue  to  vary  outside  the  boundary 
layer  in  conformity  with  the  calculated  variation  of  pressure  in 
the  inviscid  flow.  The  velocity  profiles,  and  the  profiles  of 
turbulent  kinetic-energy  shown  in  Fig.  28,  are  also  in  good 
agreement  with  the  experimental  data  except  in  the  wall  region 
near  the  tail  ( X/L  >  0.95,  say),  where  the  boundary  layer  becomes 
thick.  The  larger  values  of  mean  velocity  in  the  wall  region  of 
the  thick  boundary  layer  are  presumably  related  to  an  over¬ 
estimation  of  the  eddy-viscosity  by  the  basic  k-e  model.  The 
wall-shear  velocity  Ut  is  also  over-predicted  over  the  last  10 
percent  of  the  body  length  (Fig.  29).  In  view  of  the  observed 
disagreement,  the  measured  velocity  distributions  were  reanalyzed 
using  the  method  of  Clauser  to  determine  the  wall  shear  stress 
(see  Fig.  30).  As  shown  in  Fig.  29,  however,  the  difference 
between  these  and  the  Pres  ton- tube  values  quoted  by  the 
experimenters  is  too  sr.all  to  account  for  the  disagreement 
between  the  calculations  and  data  over  the  tail  region.  In  the 
upstream  thin  boundary-layer  region,  the  present  calculations 
also  predict  somewhat  higher  values  of  Ut.  This  may  be  due  to 
the  use  of  the  simple  1/7  -  power  law  velocity  profile  at  the 
initial  station  and  possibly  an  insufficient  grid  resolution. 

The  calculations  for  the  second  axisymmetric  body,  namely 
Afterbody  1  (L=10.06  ft),  were  performed  with  60  stations  in  the 
domain  0.50  <  X/L  <  16.25,  and  19  points  between  the  body  and  the 
external  boundary  at  R/L  =  0.8137.  The  coordinates  for  this  case 
are  shown  in  Fig.  31,  These  indicate  somewhat  milder  curvature 
changes  in  the  longitudinal  direction  compared  to  those  for  the 
Afterbody  5.  The  same  initial  conditions  as  for  Afterbody  5  were 
specified  at  the  upstream  station  X  =  0.5L,  and  a  time  increment 
t  =  1.0  was  used  again. 
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The  principal  results  of  the  calculations  for  Afterbody  1 
are  shown  in  Figs.  32  through  36.  As  can  be  seen  from  Figs.  32 
and  35,  the  solution  again  converges  in  less  than  30  global 
sweeps.  The  converged  surface  pressure  distribution  shown  in 
Fig.  32  is  in  excellent  agreement  with  the  experimental  data. 
The  pressure  along  the  wake  centerline  (X/L  >  1.0),  however, 
decays  somewhat  faster  than  the  data  in  the  near  wake  and  becomes 
slightly  negative  before  gradually  recovering  to  the  uniform 
ambient  pressure  in  the  far  wake.  The  detailed  pressure  varia¬ 
tions  in  the  radial  direction  shown  in  Fig.  33,  are  also  in  good 
agreement  with  the  data.  Note  that  the  pressure  changes  are 
relatively  mild  compared  with  the  greater  changes  observed  on 
Afterbody  5  (Fig.  26).  This  is  obviously  due  to  the  smoother 
curvature  variations  in  the  present  case.  It  is  also  interesting 
to  note  that  the  uniform  ambient  pressure  is  again  recovered  in 
the  outer  region,  and  the  viscous-inviscid  interaction  is  res¬ 
tricted  to  radial  distances  of  the  order  of  0.351.  from  the  body 
surface  and  wake  centerline. 

Detailed  comparisons  between  the  calculated  and  experimental 
mean-velocity  profiles  and  turbulent  kinetic-energy  distributions 
are  shown  in  Figs.  34  and  35,  respectively.  It  is  seen  that  the 
boundary  layer  thickness  and  half-width  of  the  wake  are  correctly 
predicted.  The  axial  (U)  and  radial  (V)  components  of  velocity 
in  the  stern  and  near  wake  region  are  also  in  very  good  agreement 
with  the  corresponding  data.  The  calculated  turbulent  kinetic- 
energy  distributions,  however,  are  somewhat  larger  in  the  thick 
boundary  layer  region.  As  expected  from  the  velocity  distribu¬ 
tions,  the  wall-shear  velocity  shown  in  Fig.  36  is  in  fairly  good 
agreement  with  the  data  even  in  the  thick  boundary  layer  region 
near  the  tail.  Here  again,  the  experimental  values  obtained  by 
Preston  tubes  agreed  reasonably  well  with  those  determined  from 
the  Clauser  plots  (Fig.  37). 

A  comparison  between  the  calculations  for  Afterbody  1  and 
Afterbody  5  indicates  that  quite  satisfactory  results  are  ob- 
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tamed  in  both  cases.  However,  some  discrepancies  are  observed 
near  the  tail  of  Afterbody  5,  which  contains  an  inflexion  point 
in  the  longitudinal  profile  and  leads  to  rather  larger  adverse 
pressure  gradients.  The  differences  in  the  boundary- layer 
behavior  in  the  tail  region  of  the  two  bodies  are  also  apparent 
from  the  Clauser  plots  shown  in  Figs.  30  and  37.  It  is  seen  from 
Fig.  37  that,  under  the  moderate  pressure  gradients  on  Afterbody 
1,  the  velocity  profiles  indicate  larger  regions  of  logarithmic 
behavior.  On  the  other  hand,  the  Clauser  plots  for  Afterbody  5 
shown  in  Fig.  30  indicate  significant  departures  from  the  law  of 
wall  and  greater  scatter.  The  basic  k-e  turbulence  model 
employed  here  does  not  include  any  corrections  for  transverse  or 
longitudinal  curvatures,  but  a  pressure-gradient  correction  has 
been  employed  in  the  wall  functions.  The  results  indicate  that 
the  latter  may  not  be  entirely  satisfactory  in  strong  pressure 
gradients.  Also,  the  somewhat  larger  turbulent  kinetic-energy 
levels  predicted  in  the  thick  boundary  layer  region  of  Afterbody 
S  suggest  that  some  mod  if i  ations  may  be  required  in  the 
turbulence  model. 

Results  similar  t  '  those  presented  here  have  also  been  ob¬ 
tained  by  a  variety  A  other  methods,  ranging  from  interacting 
bour-da  ry- laye  r  analysis  (Patel  and  Lee,  1977;  Lee,  1978;  Huang, 
Bantelli  and  Belt,  1979;  Dietz,  1980;  Hoffman,  1980),  through 
pa  r t  i  a  1  1 y-parabo  1  i  c  solutions  (Muraoka,  1980;  Hoffman,  1980; 
Hogan,  19^3),  to  a  f  u  1  ly-e  1  1  ipt.  ic  solution  (Zhou,  1982  ).  Unfor¬ 
tunately,  only  Huang  et  al.  and  Hogan  have  presented  calculations 
for  the  two  todies  cons  i  dt?red  here.  Somewhat  improved  velocity 
profiles  for  Afterbody  5  were  obtained  by  both  in  the  wall  re¬ 
gie-ins  by  incorporating  an  empirical  correction  in  the  turbulence 
model s . 

A  particularly  noteworthy  feature  of  the  present  solutions 
is  the  behavior  of  the  pressure  distribution  close  to  the  tail. 
Figures  2S  and  32  indicate  a  region  of  pressure  decrease  followed 
by  an  increase  over  the  extreme  tail  a md  then  a  monotonic  decease 


in  the  near  wake.  In  the  interacting  boundary-layer  solutions  of 
Huang  et  al.,  on  the  other  hand,  the  calculated  pressure  de¬ 
creases  continuously  at  the  tail.  We  believe  that  the  features 
predicted  by  the  present  method  are  genuine  and  are  associated 
with  the  rapid  changes  in  the  geometry  near  the  tail  as  well  as 
the  upstream  influence  of  the  complex  pressure  interaction  in  the 
tail  region.  Note  that  such  a  pressure  behavior  is  also  pre¬ 
dicted  by  the  present  and  previous  calculations  for  even  the 
simplest  case  of  the  trailing  edge  of  a  flat-plate  (Figs.  11  and 
15)  . 

Finally,  it  should  be  remarked  that  the  present  calculations 
for  Afterbody  1  and  Afterbody  5,  as  well  as  those  for  three  other 
bodies,  for  which  extensive  data  are  available  in  the  thick  boun¬ 
dary  layer,  extend  to  quite  large  distances  (~  15L)  into  the  far 
wake.  A  critical  examination  of  all  the  results  will  provide 
information  not  only  on  the  performance  of  the  turbulence  model 
in  the  presence  of  strong  lateral  and  longitudinal  surface  curva¬ 
tures  over  the  tail  region  but  also  on  the  flow  in  the  near  and 
far  axisymmetric  wakes,  for  which  the  available  experimental  data 
is  rather  limited. 

V.3  Three-Dimensional  Flow  Over  the  Stern  and  in  the  Wake  of 
Ship-Like  Bodies 

Huang,  Groves  and  Belt  (1983,  1984)  and  Groves,  Belt  and 
Huang  (1982)  have  reported  extensive  measurements  in  the  stern 
boundary  layer  of  two  ship-like  bodies  whose  cross-sections  are 
elliptic,  with  axes  ratio  3:1  and  2:1.  These  bodies  have  the 
same  lengths  and  cross-sectional  area  distributions  as  the  parent 
axisymmetric  (1:1  axes-ratio)  Afterbody  1.  We  shall  present  some 
aspects  of  the  calculations  performed  with  the  present  method  for 
these  bodies. 

The  three-dimensional  calculations  for  both  bodies  were 
performed  with  60  stations  in  the  longitudinal  direction  between 
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0.5  <  X/L  <  16.25.  Of  these,  40  are  placed  in  the  boundary-layer 
region  0.5  <■  X/L  <  1  to  capture  the  flow  features  associated  with 
the  rapid  changes  in  body  curvatures  in  the  longitudinal  and 
circumferential  directions.  The  remaining  20  stations  are  lo¬ 
cated  in  the  wake  upto  X/L  =  16.25.  Two  nonuniform  grid  nets  of 
19  x  7  and  19  x  12  mesh  points,  respectively,  are  used  at  each 
cross-section,  the  outer  boundary  being  a  cylindrical  surface  of 
constant  radius  R  =  0.8137L.  Some  views  of  the  numerically-gen¬ 
erated  coordinates  for  the  3:1  body  with  (60,  19,  7)  mesh  points 
in  the  (C,n,S)  directions  were  shown  in  Fig.  5.  Similar  coordi¬ 
nates  were  used  for  the  second  body.  For  this  case,  calculations 
were  also  performed  with  the  finer  (60,  19,  12)  mesh  to  study  the 
grid-dependence  of  the  solutions. 

The  solutions  were  started  by  specifying  the  velocity  and 
turbulence-parameter  profiles  at  the  upstream  station  where  the 
boundary  layer  is  thin  and  a  constant  ambient  pressure  was  again 
assumed  throughout  the  solution  domain.  For  simplicity,  flat- 
plate  boundary-layer  profiles  with  constant  boundary-layer  thick¬ 
ness  around  the  girth  were  specified  at  the  upstream  station, 
although  more  appropriate  initial  conditions  can  be  obtained  from 
three-dimensional  thin  boundary-layer  solutions.  The  unsteady 
marching  technique  is  again  employed  with  large  time  increments 
of  t  =  1.0  and  0.5  for  the  coarse  and  fine  mesh  calculations, 
respectively.  Plane  of  symmetry  boundary  conditions  were  spec¬ 
ified  at  9  =  0°  (?  =  2)  and  9  =  ( c  =  N)  as  follows: 

Along  9  =  0(^=2): 


r(5,n,l)  = 
9(5, n,l)  = 
4>  (  5  ,n  ,1 )  = 
W(5,n,l)  = 


r  (  5  ,  n  ,  3  ) 
0(5, n, 3) 

4>(  5  ,  n  , 3  )  ,  t 
-  W( ( 5,n  ,2) 


~  U  ,  V  ,  k  ,  e  ,  p 


and,  along  9=90° (?=N) 


r(5,n,N+l) 
9  (  S  ,  n  ,N+1 ) 
♦(«, nfN+l) 
w(  € , n ,N)  = 


r(  C , n ,N-1 ) 
it-0(  £ ,  n  ,N— 1 ) 
♦  U,n,N-l)  , 
W(C,n,N-l) 


=  U , V, k , e ,p 


where  N=6  and  11  for  the  coarse  and  fine  meshes,  respectively. 

Figures  38  and  39  show  the  convergence  histories  of  the  wall 
and  wake-centerline  pressure  distributions,  Cp  =  2p  and  the 
friction  velocity,  UT  ,  respectively,  along  6  =  90°  for  the  3:1 
elliptic  body.  It  can  be  seen  that  both  converged  monotonically 
in  about  30  time  steps  or  global  sweeps.  The  converged  pressure 
distributions  along  9=0°  and  90°,  shown  in  Fig.  40,  are  in 
fairly  good  agreement  with  the  data.  As  in  the  previous  two 
examples,  significant  pressure  variations  are  again  observed  in 
the  tail  region. 

Although  there  is  no  data  available  for  direct  comparison, 
the  pressure  distributions  across  the  boundary  layer  and  in  the 
external  flow  are  shown  in  Fig.  41  at  several  streamwise 
stations.  It  is  seen  that  the  uniform  ambient  pressure  is 
recovered  in  the  outer  region,  and  the  viscous-inviscid  interac¬ 
tion,  which  is  confined  to  radial  distances  of  the  order  of  0.4L, 
is  again  captured.  This  is  particularly  encouraging  considering 
the  practical  difficulties  of  matching  separate  viscous  and 
inviscid  solutions  for  three-dimensional  thick  boundary  layers 
and  wakes. 

Figure  42  shows  the  calculated  and  experimental  wall  pres¬ 
sure  distributions  around  the  girth  at  several  streamwise  sta¬ 
tions.  Higher  wall  pressures  are  observed  along  0  =  90°  as  a 
result  of  more  rapid  curvature  changes  along  that  line.  The 
calculations,  however,  predict  a  somewhat  stronger  three-dimen¬ 
sional  effect  as  evidenced  by  the  larger  girthwise  variation  of 
pressure.  It  is  somewhat  surprising  to  note  that  larger  discre- 
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pancies,  in  fact,  are  observed  along  0=0°  where  the  boundary 
layer  remains  relatively  thin  throughout  the  stern  region.  It  is 
also  interesting  to  note  that  the  potent ial-f low  calculat ions , 
provided  by  Huang  of  the  DTNSRDC  (David  W.  Taylor  Naval  Ship 
Research  and  Development  Center),  also  predict  lower  wall  pres¬ 
sures  along  0  =  0°,  as  shown  in  Fig.  43.  A  comparison  between 
the  potential-flow  solutions  and  the  present  calculations 
indicates  that  the  viscous-inviscid  interaction  is  confined 
largely  to  the  vicinity  of  0  =  90°  where  the  boundary  layer 
thickness  is  large.  The  large  transverse  and  longitudinal 
inviscid  pressure  gradients  around  0  =  90°  also  suggest  that 
conventional  boundary-layer  calculations  based  on  the  potential- 
flow  pressure  distribution  will  be  in  error,  and  may  predict 
premature  separation,  in  this  region. 

Figure  44  shows  the  distributions  of  the  two  nonzero 
components  (U,V)  of  mean  velocity  within  the  boundary  layer  and 
in  the  wake  in  the  two  planes  of  symmetry.  It  should  be  empha¬ 
sized  that  this  is  a  close-up  view  of  the  viscous  region  since 
the  solution  domain  extends  well  into  the  external  inviscid  flow 
field.  The  calculations  are  seen  to  be  in  good  agreement  with 
the  measurements,  which  extend  only  upto  X/L  =  0.954.  In  fact, 
the  level  of  agreement  is  comparable  to  that  achieved  in  the 
axisymmetric  flow  on  Afterbody  1  (Fig.  34).  It  is  interesting  to 
observe  the  evolution  of  the  three-dimensionality  of  the  flow 
over  the  stern.  The  boundary-layer  thickness  is  nearly  constant 
around  the  girth  at  X/L  =  0.719  but  the  thickness  along  0  =  90° 
is  almost  three  times  that  along  0  =  0°  by  X/L  =  0.954.  Although 
no  data  are  available  in  the  wake,  the  calculations  indicate  that 
the  three-dimensionality  of  the  wake  persists  for  quite  large 
distances  downstream.  However,  due  to  turbulent  diffusion  in  the 
circumferential  direction,  the  solutions  eventually  lead  to  an 
axisymmetric  far  wake.  Also,  we  note  that  conventional  boundary- 
layer  calculations  performed  by  Huang,  Groves  and  Belt  (1983)  and 
Patel,  Sarda  and  Shahshahan  (1983)  for  this  body  breakdown  around 


X/L  =  0.93  due  to  the  thickening  of  the  boundary  layer,  whereas 
the  present  method  does  not  encounter  such  a  difficulty  and 
continues  smoothly  into  the  wake. 

Figure  45  shows  the  turbulent  kinetic-energy  distributions 
along  9=0°  and  90°  at  several  streamwise  locations.  Since  not 
all  normal  Reynolds-stress  components  were  measured  on  the  planes 
of  symmetry,  the  calculated  turbulent  kinetic-energy  distribu¬ 
tions  along  9=0°  and  9  =  90°  are  compared  with  the  measurements 
at  9  =  67°  and  9  =  87°,  respectively.  It  is  seen  that  the  agree¬ 
ment  between  the  calculations  and  experiments,  in  general,  is 
quite  good  except  at  9  =  90°  and  X/L  =  0.934  ,  where  the  experi¬ 
ments  indicate  a  sudden  increase  of  turbulent  kinetic-energy, 
while  the  calculations  predict  a  continuous  decrease  of  maximum 
turbulent  kinetic-energy,  and  a  continuous  thickening  of  boundary 
layer  toward  the  tail.  The  disagreement  between  the  calculated 
and  measured  kinetic-energy  distributions  in  the  outer  part  of 
the  thick  boundary  layer  along  9  =  90°  is  also  noteworthy  since 
it  may  be  related  to  a  deficiency  of  the  turbulence  model. 

Finally,  the  wall  shear  stresses  along  9=0°  and  9  =  90°, 
which  were  not  measured  in  this  experiment,  were  determined  from 
Clauser  plots  of  the  measured  velocity  profiles  as  shown  in  Fig. 
46.  In  Fig.  47,  the  calculated  wall-shear  velocities  are 
compared  with  those  obtained  from  the  Clauser  plots.  It  is  seen 
that  the  calculations  predict  the  overall  features  observed  in 
the  experiments  but  the  values  are  somewhat  higher.  It  is  also 
interesting  to  note  that  the  present  calculations  predict  a 
higher  shear  stress  along  9  =  90°  in  the  thin  boundary  layer 
region  upto  X/L  =  0.7  even  though  a  constant  shear  stress  was 
specified  around  the  girth  at  the  upstream  station.  This 
behavior  is  also  evident  from  the  Clauser  plots  of  the  measured 
velocity  profiles  at  X/L  =  0.719. 

For  the  2:1  axes-ratio  body,  the  calculations  were  performed 
using  both  the  coarse  (60x19x7)  and  the  finer  (60x19x12)  grid. 
Figure  48  shows  the  calculated  and  experimental  pressure  distri- 


butions  on  the  wall  and  along  the  centerline  of  the  wake  on  the 
two  planes  of  symmetry.  It  is  seen  that  both  grids  give  essen¬ 
tially  the  same  results  which  are  in  good  agreement  with  the 
data.  Note  that  the  pressure  distributions  show  a  behavior  very 
similar  to  those  obtained  for  the  3:1  body  (Fig.  40)  and  the 
axisymmetric  Afterbody  1  (Fig.  32). 

The  pressure  distributions  across  the  boundary  layer  and  in 
the  external  inviscid  flow  are  compared  with  the  experimental 
data  in  Fig.  49.  The  agreement  between  the  calculations  and 
measurements,  in  general,  is  good  within  the  boundary  layer. 
However,  the  calculated  pressures  recover  slowly  to  the  uniform 
ambient  pressure  in  the  outer  region,  while  the  measurements 
indicate  an  unreal ist ical ly  abrupt  return  to  zero  pressure  just 
outside  the  boundary  layer.  This  is  particularly  evident  from 
the  close-up  views  of  the  pressure  distributions  shown  in  Fig. 
50.  The  measurements  extend  only  upto  a  normal  distance  of  the 
order  of  10  inches,  and  indicate  rather  large  normal  pressure 
gradients  at  the  outermost  measurement  station.  This  is  in  sharp 
contrast  with  the  measurements  on  the  axisymmetric  Afterbody  1 
(Fig.  33)  which  clearly  indicated  a  much  larger  region  of  normal 
pressure  variation  and  a  gradual  approach  to  the  ambient 
pressure.  In  fact,  the  calculations  for  all  three  bodies  (1:1, 
2:1  and  3:1)  indicate  about  the  same  extent  of  viscous- inviscid 
interaction,  and  therefore  the  disagreement  seen  in  Figs.  49  and 
50  outside  the  boundary  layer  is  most  likely  due  to  systematic 
errors  in  the  measurement  of  the  rather  small  pressure  changes. 
It  should  be  remarked  here  also  that  the  measured  wall  pressures 
appear  to  be  inconsistent  with  the  pressures  measured  in  the  flow 
field.  In  particular,  the  wall  pressures  at  X/L  =  0.858  and 
0.894  along  9  =  90°  are  considerably  lower  than  those  measured 
very  close  to  the  wall  (see  Fig.  50b).  It  is  also  of  interest  to 
note  that  the  present  calculations  predict  a  decrease  of  pressure 
in  the  near-wall  region  along  0=0°  upto  X/L  =  0.95,  whereas  the 
experiments  indicate  a  continuous  increase  toward  the  body 


98 


surface  in  the  tail  region.  If  these  differences  are  taken  into 
account,  the  general  level  of  agreement  between  the  calculations 
and  experiments  is  actually  much  better  than  that  indicated  by 
the  wall  pressure  distribution  shown  in  Fig.  48. 

For  further  understanding  of  the  three-dimensional  pressure 
field,  the  calculated  and  experimental  wall  pressures  around  the 
girth  are  shown  in  Fig.  51.  The  calculations  are  in  fairly  good 
agreement  with  data  and  show  a  behavior  very  similar  to  that 
observed  on  the  3:1  body,  except  that  the  girthwise  pressure 
changes  are  somewhat  smaller,  as  expected.  It  should  be  re¬ 
marked,  however,  that  considerable  scattering  of  the  data  is  seen 
near  0  =  90°  while  the  calculations  predict  a  smooth  variation 
toward  the  symmetry  plane.  Nevertheless,  the  overall  agreement 
in  the  pressure  distributions  is  quite  satisfactory. 

Detailed  comparisons  of  the  calculated  velocity  profiles 
with  the  experimental  data  are  shown  in  Figs.  52  and  53  for  the 
coarse  and  finer  grid  calculations,  respectively.  Both  calcula¬ 
tions  are  in  good  agreement  with  the  corresponding  data  in  the 
boundary-layer  region  (X/L  <  1),  but  the  calculations  with  the 
finer  gird  predict  a  somewhat  stronger  three-dimensionality  in 
the  far  wake,  presumably  due  to  the  better  resolution  of 
travsverse  gradients  near  0  =  0°.  The  evolution  of  the  three- 
dimensionality  of  the  flow  over  the  stern  is  again  well 
predicted. 

Figure  54  shows  a  comparison  between  the  calculated  wall- 
shear  velocities  along  0=0°  and  90°  and  those  measured  by  Huang 
et  al.  using  Preston  tubes.  It  is  seen  that  the  calculations 
exhibit  the  same  trends  as  the  experimental  data,  but  the 
calculated  values  are  much  higher.  This  is  particularly 
surprising  in  view  of  the  fact  that  the  calculated  velocity  pro¬ 
files  are  in  fairly  good  agreement  with  the  data.  In  order  to 
resolve  this  apparent  discrepancy,  the  experimental  velocity 
profiles  were  again  analyzed  using  Clauser  plots,  which  are  shown 
in  Fig.  55.  The  wall-shear  velocities  determined  in  this  manner 


are  shown  in  Fig.  56.  These  are  in  much  better  agreement  with 
the  calculations. 

In  Fig.  57,  the  calculated  turbulent  kinetic-energy  distri¬ 
butions  are  compared  with  the  experimental  data  at  several 
streamwise  stations.  In  the  thick  boundary  layer  region,  the 
calculated  energy  levels  are  comparable  with  those  obtained  for 
the  3:1  body,  but  are  much  higher  than  the  data  on  the  2:1 

body.  This  is  somewhat  surprising  in  view  of  the  good  agreement 
observed  in  the  previous  case  (see  Fig.  45). 

It  is  instructive  to  compare  the  calculations  and  measure¬ 
ments  at  a  fixed  streamwise  station,  X/L  =  0.934  say,  common  to 

all  three  bodies  of  the  same  family.  in  the  axisymmetric  case. 

Fig.  35  shows  that  the  maximum  measured  value  of  k  is  of  the 
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order  of  2.9  x  10  and  this  compares  favorably  with  the 

calculated  values.  On  the  2:1  and  3:1  bodies  (Figs.  57  and  45, 

respectively)  the  maximum  measured  values  in  the  relatively  thin 
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boundary  layer  along  the  flatter  surface  at  9  =  0°  are  3.4  x  10 

__  3 

and  4.0  x  10  ,  and  the  corresponding  calculations  indicate  3.9  x 

10  3.  This  level  of  agreement  is  considered  quite  satisfactory 
and  the  general  increase  in  k  with  increasing  axes  ratio  observed 
in  the  experiments  and  predicted  by  the  calculations  may  be  as¬ 
sociated  with  the  divergence  of  the  mean  streamlines  out  of 
the  0=0°  plane  of  symmetry.  The  situation  along  the  9  =  90° 
plane,  where  the  boundary  layer  is  thick  and  the  transverse  cur¬ 
vature  is  large,  is  more  complex.  The  measured  maximum  values  of 
k  are  1.5  x  10  ^  for  the  2:1  body  and  3.5  x  10-^  for  the  3:1  body 
while  the  calculations  indicate  2.4  x  10-^  and  2.0  x  10~^, 
respectively,  i.e.,  a  systematic  decrease  from  the  axisymmetric 
case.  As  noted  earlier,  the  high  measured  value  at  X/L  «  0.934 
on  the  3:1  body  (see  Fig.  45)  is  not  consistent  with  the 
measurements  further  upstream.  If  this  is  disregarded,  the 
measurements  and  calculations  depict  similar  trends  with  respect 
to  the  maximum  turbulence  levels.  However,  the  rather  dramatic 
difference  in  the  detailed  shapes  of  the  measured  k  profiles  on 
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the  2:1  and  3:1  bodies  (compare  Figs.  57  and  45)  along  0  = 
remains  a  mystery.  The  calculations  do  not  indicate  any  major 
changes  in  the  profile  shape  among  the  three  bodies  although,  in 
all  three  cases,  the  predictions  indicate  greater  diffusion  in 
the  outer  parts  of  the  boundary  layer.  The  latter  feature  ap¬ 
pears  to  be  a  deficiency  of  the  turbulence  model.  Obviously,  a 
much  more  detailed  analysis  of  the  experimental  and  calculated 
results  than  attempted  here  is  required  in  order  to  gain  a  better 
insight  into  the  performance  of  the  turbulence  model. 

Thus  far,  we  have  considered  the  details  of  the  flow  in  the 
planes  of  symmetry  of  the  three-dimensional  bodies  since  direct 
comparison  is  possible  without  any  interpolation  in  the  data  or 
calculations.  Comparisons  off  the  symmetry  planes  are  not  as 
straightforward  due  to  the  quite  different  coordinates  utilized 
in  the  experiments  and  calculations.  An  indication  of  the  over¬ 
all  three-dimensionality  of  the  flow  on  the  2:1  body  is,  however, 
provided  by  plotting  the  velocity  vectors  in  several  cross-sec¬ 
tions.  Since,  in  the  present  formulation,  the  longitudinal  com¬ 
ponent  (U)  is  parallel  to  the  body  axis,  it  is  necessary  to  show 
only  the  resultant  of  the  radial  (V)  and  circumferential  (W) 
components.  The  results  of  the  finer-grid  solutions  are  shown  in 
this  manner  in  Fig.  58.  Note  that  these  plots  depict  only  a 
portion  of  the  solution  in  the  vicinity  of  the  body,  which  is 
also  shown  to  scale. 

In  the  thin  boundary  layer  region  upto  X/L  =  0.853  (Fig. 
58a, b),  the  velocity  vectors  are  essentially  normal  to  the  wall, 
and  directed  towards  the  wall  due  to  the  longitudinal  curvature 
of  the  surface.  The  crossflow  angles  are  relatively  small.  The 
increasing  normal  velocity  components  predicted  in  the  thick 
boundary  layer  region  near  the  tail  (Figs.  58c, d,e)  are  obviously 
due  to  the  rapid  changes  of  body  shape.  Although  the  magnitude 
of  the  secondary  flow  is  quite  small,  the  flow  vectors  indicate 
the  development  of  a  small  zone  of  circumferential  flow  reversal 
or  vortical  flow  in  the  tail  region.  This  is  particularly  evi- 
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dent  in  the  wake  region  where  the  radial  velocity  components 
become  small  also.  For  a  better  understanding  of  the  development 
of  this  feature,  the  inner  regions  of  the  flow  at  a  few  stations 
near  the  tail  of  the  body  are  further  enlarged  and  shown  in  Fig. 
59.  A  longitudinal  vortex  like  structure  is  clearly  evident  in 
the  near  wake  (Fig.  59c).  Its  size  appears  to  grow  but  the 
strength  diminishes  with  downstream  distance.  The  results  shown 
in  Figs.  58  and  59  indicate  that  the  flow  on  the  2:1  elliptic 
body  is  only  weakly  three  dimensional.  A  similar  conclusion  can 
be  drawn  also  for  the  3:1  body. 

Finally,  it  should  be  emphasized  that  the  present  calcula¬ 
tions  utilize  the  radial  (V)  and  circumferential  (W)  components 
of  velocity  in  a  basic  cylindrical  polar  coordinate  system  re¬ 
ferenced  to  the  body  axis,  and  not  the  conventional  boundary- 
layer  velocity  components  (i.e.,  normal  and  tangential  to  the 
body  surface).  Therefore,  the  calculated  V  and  W  component 
velocities  are,  in  fact,  quite  large  and  of  the  same  order  of 
magnitude  even  though  the  magnitude  of  the  secondary  flow  in 
boundary-layer  coordinates  is  relatively  small.  The  calculations 
described  here  presented  no  difficulty  in  handling  these  large 
component  velocities,  and  the  numerical  method  is  not  restricted 
to  small  crossflows. 


VI.  CONCLUSIONS 


A  versatile  new  method  for  the  calculation  of  the  flow  over 
the  stern  and  in  the  wake  of  ship-like  bodies  has  been  described 
and  several  examples  have  been  presented  to  demonstrate  its  capa¬ 
bilities.  The  method  solves  the  partially-parabol ic  (or  "semi- 
elliptic"  or  "parabolized")  Reynolds-averaged  Navier-Stokes  equa¬ 
tions.  Among  the  distinctive  features  of  the  method  are  the 
following : 

(1)  Numerically-generated  nonorthogonal  coordinates  are  used 
to  facilitate  application  of  the  method  to  a  wide 
variety  of  body  shapes. 

(2)  The  finite-analytic  scheme,  which  uses  analytic 

solutions  of  the  locally-linearized  convective  transport 
equations,  enables  the  resolution  of  the  flow  in  regions 
of  larger  velocity  gradients  with  comparatively  few  grid 
points  in  the  solution  domain  and  results  in  an 
economical  method. 

(3)  The  transport  equations  are  solved  by  a  space-  and  time¬ 
marching  technique,  which  enables  future  extensions  to 
unsteady  and  fully-elliptic  problems. 

(4)  The  elliptic  pressure  field  is  determined  by  a  new  two- 
step  global  iteration  technique  which  results  in  rapid 
and  monotonic  convergence  of  the  entire  solution.  This 
algorithm  also  facilitates  the  future  extension  of  the 
method  to  fully-elliptic  equations. 

(5)  The  method  enables  the  use  of  a  large  solution  domain 
and  requires  only  the  body  geometry  and  initial 
boundary-layer  information  as  inputs.  Thus,  the  need 
for  separate  viscous  and  inviscid  solutions,  and 


matching  between  them,  is  eliminated. 

(6)  For  turbulent  flows,  the  standard  k-e  model  has  been 
used.  However,  the  usual  sensitivity  of  the  solutions 
to  the  near-wall  grid  distribution  has  been  greatly 
reduced  by  anchoring  the  solutions  to  the  wall  functions 
at  two  grid  nodes  adjacent  to  a  solid  wall. 

All  calculations  presented  here  were  performed  on  a  Prime 
750  minicomputer  with  computing  times  of  less  than  10  minutes  for 
the  two-dimensional  and  axisymmetric  cases,  and  of  the  order  of 
two  hours  for  three-dimensional  bodies.  It  should,  therefore,  be 
possible  to  use  the  method  for  practical  appl icat ions . 

While  the  overall  numerical  framework  of  the  method  has  been 
established  and  tested,  and  calculations  for  practical  hull  forms 
are  now  in  progress,  there  are  several  basic  physical  and  nu¬ 
merical  aspects  of  the  method  which  require  further  scrutiny  and 
development.  Among  these  are  the  following: 

(1)  The  influence  of  the  grid  distribution  on  the  solution 
of  the  flow  equations:  more  general  grid-control  func¬ 
tions  need  to  be  identified  to  provide  the  greater  flex¬ 
ibility  required  for  complex  ship  geometries. 

(2)  The  influence  of  the  initial  conditions:  detailed  three- 
dimensional  boundary-layer  calculations  may  be  required 
at  the  upstream  station  to  provide  more  accurate  initial 
conditions  for  complex  three-dimensional  body  geome- 


(3)  The  influence  of  the  choice  of  coordinates  on  the  solu¬ 


tions:  For  some  applications,  and  extensions  of  the 
method  to  treat  separated  flows,  it  may  be  necessary  to 
use  velocity  components  aligned  with  the  coordinate 

lines.  The  necessary  equations  have  been  derived. 
These  contain  many  more  geometrical  parameters  which 
need  to  be  either  stored  or  recalculated. 

(4)  The  sensitivity  of  the  solutions  to  the  turbulence 

model:  the  validity  of  the  basic  k-e  turbulence  model 
employed  here,  in  general,  and  the  limitations  of  the 
wall  functions,  in  particular,  in  complex  three-dimen¬ 
sional  flows  needs  to  be  studied  by  more  critical  ex¬ 
amination  of  the  available  experimental  data  and  corre¬ 
sponding  calculations. 

(5)  The  limitations  of  the  part ial ly-parabol ic  approxima¬ 

tions:  this  can  be  investigated  by  comparison  with  the 
corresponding  f u 1 ly-el 1 ipt ic  solutions. 

Finally,  we  note  that  the  generality  of  the  approach  pre¬ 
sented  here  offers  encouraging  prospects  for  incorporating  pro- 
pellors  and  appendages  in  the  solution  domain  in  a  more  compre¬ 
hensive  manner  than  has  been  possible  with  presently  available 
methods.  Also,  the  structure  of  the  solution  algorithms  selected 
is  such  that  the  method  can  be  extended  to  treat  ful ly-el 1 ipt ic , 
unsteady,  three-dimensional  flows. 
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APPENDIX  I.  COORDINATE  TRANSFORMATIONS 

This  appendix  summarizes  the  relationships  between  a  basic 
orthogonal  coordinate  system  x1  (x'*',  x^,  x^)  and  a  general  non- 
orthogonal  system  ^ )  ,  These  relations  are  used  in  the 

text  to  invert  the  Poisson  equations  for  numerical  grid  genera¬ 
tion  and  to  transform  the  Reynolds  equations  from  cylindrical 
polar  coordinates  (x,  r,  0)  to  the  nonorthogonal  numerical 
coordinates  (£,  n,  c). 

In  the  general  coordinates  £ 1  ,  any  line  element  vector  dr  is 
related  to  the  differential  change  of  coordinate  dS1  by 


3r 

/v 

dr  =  - r  d5 

~  1 


( A—  1  ) 


where  r  is  the  position  vector,  and  the  expression  is  summed  over 

t^e  three  values  of  i.  It  is  convenient  to  define  the  derivative 

- -  as  a  covariant  base  vector  a.,  i.e., 

, .  i  ~i 


a  .  = 

~i 


( A- 2 ) 


so  that  the  line  element  vector  can  be  written 


dr  =  dS  a. 

~  ~i 


(A-3) 


More  generally,  any  vector  A  other  than  a  position  vector  r  can 
be  decomposed  in  the  same  manner  as 


&  =  A  a. 


( A- 4  ) 


where  A1  denote  the  contravar iant  components  of  vector  A  . 

The  reciprocal  ( contravar iant)  base  vectors  are  defined  by 


.  a  .  x  a, 
l  ~  i  ~k 

3  =  L~i - 


_  1  m  n 

=  —  c  „  a  .  a. 

J  Zmn  j  k 


( A— 5 ) 


with  i,  j,  k  in  cyclic  order,  so  that 


AI-1 


( A- 6 ) 


1  1  t  1 

a  *  3j  =  3j  •  3  =  6j 

where  the  Kronecker  delta  61  is  a  mixed  second-rank  tensor  de- 

J 

fined  by 


i  \ lf  1  =  j 
6  1  = 

3  <  0,  i  *  j 

and  J  is  the  Jacobian  defined  by 


( A-  7) 


J  =  a. 

~  l 


( a  .  *  a.  ) 
~j  ~k 


( A-8  ) 


Any  vector  A  can  also  be  decomposed  in  terms  of  the  contra- 
variant  base  vectors  g1  as 


A  =  A.  a' 


( A-9 ) 


where  A^  are  the  covariant  components  of  vector  A.  The  covar¬ 
iant  and  contravariant  base  vectors  are  not  usually  of  unit 
length,  but  a  set  of  unit  vectors  tangent  and  normal  to  the  co¬ 
ordinate  lines  can  be  defined  as 


T  .  = 

~  1 


a  . 

~i 


a  . 


(A-lOa) 


i 

D  = 


(A-lOb) 


respectively.  The  distance  ds  between  any  two  points 

r  and  r  +  dr  is  given  by 


(ds)2  =  dr  •  dr  =  (a{  •  )  d^dC^ 


( A- 1  1  ) 


where 


s  g^dc'dC3 


AI-2 


(  A- 1  2  ) 


is  the  "metric  tensor".  It  is  clear  from  this  definition  that 
g  —  =  gj^  in  any  coordinate  system.  A  conjugate  ( contravariant ) 
metric  tensor  g1-^  can  be  defined  in  a  similar  manner, 


a' 


( A-l 3 ) 


and  is  related  to  the  metric  tensor  g—  through  equation  (A-6)  by 

g  gkj  ‘  gjkg  ■  sj  (A‘14) 

In  other  words,  the  components  of  the  conjugate  metric  tensor  g1^ 
are  the  inverse  of  those  of  Using  these  equations,  a  covar¬ 

iant  component  A^  can  be  converted  into  a  contravariant  component 
A1  by  the  relation 

At  =  gtj  A^  (A-15) 

and  vice  versa,  i.e., 

A1  =  A j  (A- 16 ) 

Also,  the  contravariant  base  vectors  a1  are  related  to  the  covar¬ 
iant  base  vectors  a.  by 

=  9 i j  a1  ( A- 1 7 ) 

and 

a1  =  g1J  a-j  <a-18) 

Finally,  we  note  that  neither  A1  nor  Aj  represent  the  physical 
components  of  the  vector  A  .  If  the  physical  components  are  de¬ 
noted  by  A(i),  they  are  given  by 


A(i)  =  (g^)3^2  A*  =  (g^)3^2  g^  A ^  (no  sum  on  i)  (A-19) 

Now,  in  the  orthogonal  coordinates  x1,  any  vector  A  can  he 
expressed  as 


A  = 


( A- 20 ) 


where  E1  and  E^  are,  respectively,  the  contravar iant  and  covar¬ 
iant  components  of  A  ,  e.  are  the  covariant  base  vectors  given  by 


( hi  i ' 


Vi- 


Vi> 


(  A-  2 1 ) 


and 


i 

e 


are 


the  contravar iant  (conjugate)  base  vectors  defined  by 


such  that 


( A- 2  2  ) 


i 


where  h^  ,  h2  and 
factors  given  by 


h^  are  the 


( A- 2  3 ) 

usual  metric  coefficients  or  scale 


( ds ) 2=  (h1dx1)2+  (h2dx2)2+  (h3dx3)2  (A-24) 

=  h .  .  dx ^dx^ 
ij 

The  metric  tensor  h—  can  also  be  obtained  from  the  following 
relation 


h  .  . 

T  T 


e  .  = 


hj  6*5*=  h26i 


3 

l 


(  A-  2  5  ) 


E 


E 


(  h  .  .  )  = 
ij 


o  h; 


o  o 


( A-25a) 


Similarly,  the  conjugate  metric  tensor  h1^  in  the  orthogonal  co¬ 
ordinates  is 


j  ..... 
.  1  j  l  1  r  1  ritl  1  a 

h  =e  •  e  -  £  — —  o  „  6-i  =  — —  6. 


1=1  h2  1  1  h2  j 


( A- 2 6  ) 


(h1^)  = 


( A-26a ) 


As  before,  the  physical  components  E(i)  of  any  vector  A  in  the 
orthogonal  coordinates  are  given  by 


E (  i )  =  I-4E1  = 


-j^—  E ^  (no  sum  on  i) 
i 


( A-27) 


Since  the  equations  of  fluid  flow  are  usually  written  in 
orthogonal  coordinate  systems  x1,  it  is  desirable  to  establish 
the  transformations  between  the  basic  x1  system  and  the  general 
nonorthogonal  system  5 1 .  From  equations  (A-3)  and  (A-20),  the 
line  element  vector  dr  in  the  two  systems  is  given  by 


dr  =  dx^e.  =  dS 1  a  . 


( A-28 ) 


S  ince 


dxj  =  lx_  dCi 


(  A-29  ) 


the  covariant  vector  a.  in  £  coordinates  is  related  to  e. 

~i 

(equation  A-21)  by 

.  i  1  't ..  2  ^ 1 


3x  ^  3x  .  3x  ,  3  x  ' 

:£•;—(  h  i  ■  ,h~  .  ,  h  ~  ■) 

3£x  J  1  3C  3£  3£ 


(A-30) 


th  l 

If  the  4  component  of  is  denoted  by  a.  ,  a  second  rank  mixed 


tensor  can  be  defined  as 

„  „  .  * 
„  Jl  _  ,  J  .  K  3x 

ai  "  (^i}  "  hi  ~~T 


( A-31  ) 
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(  A- 3  2  ) 


The  components  of  the  reciprocal  ( contravariant)  base  vector  a 
are  given  by 


i  _  / JL_  All  _L_  Ail  1_  lllx 

a  V  -  1  '  h.  .  2  '  h.  a  31 

1  3x  2  3x  3  3x 

or  by  equation  (A-5),  i.e. 

i  l  .  „  .  1  „  m  n 

3  j  (^j  =  J  *mn  aj  ak 


( A- 3  3 ) 


with  (i,  j,  k)  in  cyclic  order. 

For  convenience,  we  shall  define 


.  i  T  i 

b  =  a  .  x  a,  =  J  a 
~ J  ~k 


( A-34 ) 


so  that  b*  ,  the  ith  component  of  J  g1,  using  equation  (A-30),  is 


bi  -  <»j  x  *k>t  * 


,  ni  ,  n 

_m  _n  _  /  3x  3x 

e  „  a  •  a.  —  h  n  ( r  — — 

imn  j  k  m  n  „,k 


35J  3£ 


3xn  3xm 
3^  3Ck 


and,  from  equations  (A-6)  and  (A-34) 


l  .  i  _  ,i 
aj  b»  ■ J  5j 


(  A-  3  9  ) 


where  (l,  m,  n)  are  in  cyclic  order  also. 

The  relationships  between  the  metric  tensors  in  the  E.1  and 
x1  coordinates  can  now  be  obtained  as  follows 


g  .  .  =  a  .  •  a  .  =  £ 

i J  ~i  ~1 


2=1 


i  i  ; 
a  .  a  .  =  E 
i  1 


1  =  1 


x2  _3x^ 

1  a?1  3Sj 


( A- 4  0  ) 


The  conjugate  (contravar iant )  metric  tensor  g1^  is 

3 


l  j  l 

g  =  a 


i  _  1_ 


3 

E 


_  -  b*bj  =  -I  b*  hJ0 
j2  l=1  t  Z  g  l=l  l  l 


( A- 4 1  ) 


2 

where  g  =  J  is  the  determinant  of  the  metric  tensor  g^j* 

An  alternative  expression  for  g1J  can  be  obtained  from  equa¬ 
tion  (A-14)  by  recognizing  that  g1^  is  the  inverse  of  gjjr  i.e., 


=  'Stj1'1  *  i  <9mj  9nk  -  9mk  9nj > 


( A-4  2 ) 


with  both  (i,  j,  k)  and  (l,  m,  n)  in  cyclic  order. 

It  should  be  remarked  here  that  the  metric  tensor  g—  and 
its  conjugate  g1-^  follow  the  tensor  transformation  rules 


.  t  m 

g.  .  -  h, 

ij  em  Hi  35j 


( A-43a) 


and 


i  j  _  .  2m  3C1  3S-1 

9  -  h  7171; 

3x  3x 
2  m 


( A-43b) 


Where  h^m  and  h  are  the  metric  coefficients  defined  in  equa¬ 
tions  ( A-24 )  and  (A-25)  for  the  orthogonal  coordinates  x1.  The 


AI-7 


last  two  equations  are  identical  with  equations  (A-40)  and  (A- 
41 ) . 


In  addition  to  the  above  geometrical  relationships,  we  need 
relationships  to  transform  the  Reynolds  and  continuity  equations 
from  the  x 1  coordinates  to  the  S1  coordinates.  It  is  therefore 
necessary  to  obtain  expressions  for  the  gradient,  divergence  and 
Laplacian  operators  in  the  general  coordinates  51.  These  opera¬ 
tors  involve  the  derivatives  of  either  scalar  or  vector  quanti¬ 
ties  with  respect  to  C1.  In  Cartesian  coordinates,  the  partial 
derivatives  always  give  higher  order  tensors.  This  is,  however, 
not  true  for  general  curvilinear  coordinates.  For  example,  al¬ 
though  3 <(> / 3 ^ 1  is  a  covariant  vector  for  any  scalar  function  , 
the  second  derivative  32$/3£*3£^  does  not  give  a  covariant 
second-order  tensor.  In  this  respect,  the  curvilinear  system  is 
not  as  convenient  as  the  Cartesian  system.  It  is  important  to 
define  a  derivative  which  preserves  tensor  character  under  co¬ 
ordinate  transformations  so  that  the  effect  of  coordinate  trans¬ 
formations  on  derivatives  need  not  be  considered  explicitly.  A 
derivative  with  the  desired  tensor  transformation  properties  can 
be  defined  as  follows: 
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Thisrt  is 
A*  2*’ 

6iV- 


called  the  "covariant  derivative"  of  the  general  tensor 
’ a  i 

gr  ,  and  r  .  k  is  the  Christoffel  symbol  defined  by 
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It  is  clear  that  the  Christoffel  symbol  is  symmetric  in  j  and 
k.  Furthermore,  r1.  is  related  to  the  Jacobian  J  or  / g  by 


r1  =  r i  =  I  3  log  g  =  L_  3g 

ij  ji  2  3^i  2g  g^i 
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1 

J 


3J 
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( A— 46 ) 


since  g  =  J^. 

.  a  a  ...  a 

It  is  noted  that  both  -2—-  A«V  «r  and  r^.  are  not  ten- 
sors;  however,  the  covariant ^derivative  (A-44)  does  preserve  the 
tensor  transformation  rule  under  coordinate  transformation  and  is 
therefore  a  higher  order  tensor.  The  result  of  equation  (A-44) 
is  that  we  need  not  consider  explicitly  the  effect  of  coordinate 

transformations  on  tensor  derivatives.  For  example,  since  g.., 

i  i  ,  ^ 

and  g  J ,  are  zero  in  Cartesian  system,  they  are  also  zero  in  all 

coordinate  systems. 


From  equation  (A-44),  the  covariant  derivative  of  a  covar¬ 
iant  vector  A ^  with  respect  to  51  is  given  by 


A, 


if  1 


3A.  u 

— *  -  rk-  a 
1J  k 


(A-47) 


Also,  the  covariant  derivative  of  a  contravar iant  vector  A1 
is  given  by 


A 
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» 3 


( A-48  ) 


The  divergence  of  a  vector  A  is  a  scalar  quantity  and  is 
invariant  under  coordinate  transformation.  It  can  therefore  be 
obtained  by  a  covariant  derivative  of  the  contravar ient  component 
A1  and  can  be  written  in  any  general  non-orthogonal  coordinates 
as 


d  i  v  A 


+  r: . 
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( A-49  ) 
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Onlv  in  Cartesian  coordinates  div  A  reduces  to  - r  because  the 

Christoffel  symbol  r1,  is  then  zero.  In  other  coordinates,  the 
divergence  operator  is  given  by  eq  (A-46)  as 


d  iv  A 


( JA1 ) 


( A-50a) 


where  A1  is  related  to  the  physical  components  E (t)  in  the  or¬ 
thogonal  coordinates  x1  by  the  tensor  transformation  rule 


A 


i 


E(2  ) 


Thus 


div  A  =  ~  -K  (  b  J  EU))  ( A-50b ) 

J  3  C  1 

In  the  basic  orthogonal  coordinate  system  x1,  the  ith  com¬ 
ponent  of  the  gradient  operator  is 


1  3<t> 

(,*>i  =  tt  rr 

l  3x 

If  we  apply  the  chain  rule  for  differentiation,  and  use  equations 
( A- 33  and  34),  the  gradient  operator  in  the  general  nonorthogonal 
coordinates  can  be  written 

( A-52 ) 

scalar  function  <f>  is  defined 

v2<p  =  div  (V<t>)  ( A-53  ) 


(  V4C 
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_M  =  I  bl  3i{> . 

3x*  3^  3  i 


Finally,  the  Laplacian  V  $  of  a 


by 


3  <t>  3  9 

where  V  $  =  - r  a.  and  — 

di{  i 

coordinates  Z  .  Since 
only  to  the  contravar iant 
component  of  7$  should  be 


is  the  ith  component  of  V 4>  in  general 
e  divergence  expression  (A-50)  applies 
component,  the  associated  contravar iant 
used ,  i  .  e  .  , 
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V2*  =  div  (V*)1  =  div  (glj  -^-r)  =  ^  ^~T  (J  ^T)  (A-55a) 

a?  J  ^  *  r  1  3  r  J 


V2*  =  glj 


•  +  \  (a  gij>  ^ 


3C13CJ  J  U1 


(A-55b) 


The  general  transformation  formulas  listed  above  have  been 
used  in  the  text  to  (a)  generate  the  nonorthogonal  numerical  grid 
in  the  physical  domain  and  (b)  transform  the  governing  equations 
of  fluid  motion  from  a  basic  orthogonal  coordinate  system  into 
the  nonorthogonal  numerical  coordinates.  The  procedures  used  to 
accomplish  this  are  outlined  below. 

( a )  Numerical  grid  generation 

The  numerical  coordinates  i1  are  related  to  the  basic  or¬ 
thogonal  coordinates  x1  by  the  set  of  Poisson  equations  (equation 
III-19  in  the  text) 


vV  .  f* 


(A-56) 


4  ,  j  1 

where  V  is  the  Laplacian  in  the  x  coordinates  and  f  are  func¬ 
tions  of  C1  which  control  the  grid  spacing.  Using  equation  (A- 
55b),  this  may  be  written 


f‘  =  gU  *  A  -t-  ,J  gU, 


13  ,  _  H 

J  77  i  (J  g 


( A-57) 


Since  J  =  h^ 


i  i  Z  £ 

and  g  =  h  in  the  orthogonal  system,  f  may 


be  expressed  in  terms  of  the  metric  coefficients  as 
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51  h,h'h,  T1  t  ‘J  3>  (no  sum  on  <A'5S) 

1  2  3  3x  hj 

using  equation  (A-57),  the  second  term  on  the  right  hand  side  of 

l 

equation  (A-55)  can  be  expressed  in  terms  of  f  and  the  expres¬ 
sion  for  the  Laplacian  in  the  numerical  coordinates  governed  by 
equation  (A-56)  becomes 
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(  A- 5  9  ) 


Thus,  the  inverse  form  of  equations  (A-56)  which  can  be  used  for 
grid  generation  becomes 


„2  i 
7  x 


i  J 


a2  i 
3  x 

3  e. 1 3  4  J 


+  f 


1  Jx. 
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( A- 60 ) 


l 

where  g1-1  is  evaluated  in  terms  of  g^j  by  equation  (A-42)  and 
may  be  assigned  appropriate  values  to  control  the  grid  spacing. 
Equations  (A-60)  with  i  =  1,  2,  3  are  equations  (III-19)  in  the 
text. 


( b )  Transformation  of  the  flow  equations 

Consider  the  equations  of  motion  in  orthogonal  curvilinear 

12  3  , 

coordinates  (x  ,  x  ,  x  )  for  unsteady,  three-dimensional,  incom¬ 
pressible  flows.  The  exact  Reynolds-averaged  equations  of  con¬ 
tinuity  and  momentum  of  the  mean  flow  (Nash  and  Patel  (1972))  in 
dimensionless  form  are. 
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with 
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nr-.-  the  curvature  parameters,  a  ,  a  etc.  are  functions  only 
of  the  curvature  parameters  k—  and  their  derivatives,  and  can  be 
wr  i  t  ten 


all  = 

-  'k2i4 

9 

k12  + 

k 1 3  *  k31> 

X  - 

2  2 

-  -,2 

"  k223 

4  k  2^  4  k  1  2 

1  3  3  = 

-  ,kH 

f  k  3  l 

4  k12  4  k  2  1 

A  r  -  1  4 


'12 


a 


13 


a 


21 


a 


23 


a 


31 


a 


32 


1 

3k12 

1 

3k21 

hl 

3x 1 

‘  h2 

3x  2 

1 

3k  1 3 

1 

3k  3 1 

hl 

3X1 

"  h3 

3x3 

1 

3k21 

1 

3  k  1 2 

h2 

2  2 

3x 

"  hl 

^  1 

3  x 

1 

23 

1 

3k 

32 

h2 

3x2 

'  h3 

3x3 

1 

3k31 

1 

3k13 

h3 

3x3 

"  hl 

3X1 

1 

3k32 

1_ 

3k23 

h3 

3x3 

h2 

3x2 

k32(k31+  k21}  +  k  3 1 k  1 2 


k23(k31+  k21) 


+  k21k13 


k31(k12+  k32)  +  k  32  k  2 1 


k 1 3 ( k 1 2  +  k32}  +  k12k23 
k21  (k23+  k13)  +  k  23  k  3 1 


k 1 2 ( k  23  +  k13) 


+  k 1 3  k  32 


The  above  equations  can  be  transformed  to  the  general  cur¬ 
vilinear  coordinates  C1  by  using  the  expressions  for  divergence, 
gradient  and  Laplacian.  First,  the  continuity  equation  (A-61)  is 
transformed  from  the  coordinates  to  the  £ 1  coordinates  by 
equation  (A-50) 


div  y  =  \  - r-  (  b ,  U(  A  )  )  =  0  (A- 66) 

J  a?i 

where  U(«)  are  the  physical  components  (U,  V,  W)  of  velocity 
vector  U  in  the  orthogonal  coordinates  x1.  Equation  (A-66)  when 
written  in  terms  of  the  velocity  components  (U,  V,  W)  in  the 
cylindrical  polar  coordinates  (x,  r,  9)  is  equation  ( T 1 1  —  1 )  in 
the  text. 

Similarly,  the  momentum  equations  (A-62)  through  (A-64)  are 
transformed  to  S1  (  £  ,  n  ,  t )  coordinates  by  equations  (A-52)  and 
( A- 59  )  , 
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In  the  present  study,  the  two  equation  (k  -  e)  turbulence 

model 

together  with  the  isotropic  eddy-viscosity  assumption  is 

i 

used  to  model  the  Reynolds  stresses.  In  other  words. 

the  Rey- 

'# 

nolds 

stresses  are  assumed  to  be  related  to  the  mean 

rates  of 

strain 

by  means  of  the  turbulent  viscosity  as  follows: 
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The  eddy  viscosity  v  is  related  to  the  dimensionless  turbu¬ 
lent  kinetic  energy  k  and  the  rate  of  dissipation  of  turbulent 
energy  e  by 


v  =  c  — — 

t  p  e 


(A- 71) 


and  k  and  e  are  governed  by  convective  transport  equations  of  the 
following  form 
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where  G  is  the  turbulent  generation  term 
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and  the  effective  Reynolds  number  is  defined  as 
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Using  equation  (A- 70),  (A-71)  and  (A-75),  the  momentum  equa¬ 
tions  ( A- 6  7 )  through  (A-69)  and  the  turbulence-model  equations 
( A- 72 )  and  (A-73)  are  transformed  to 
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U  +  a . v  +  a 


+  <k  V  -  k0oW) W  +  (k 


U  U  + 


U  +  2k32W)  -  h  3  (h  2 
3  3x  2  3x 


_  2k  JL  _iw  +  ! _ 3w 

/K 32  h-  a  3  23  h0  ,  2 

3  3x  2  3x 


1  h2  3x2  +  a21U  +  a22V  +  a23W} 


3  v  3  v 
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h_  _  1  ;  h.  .  1  .  2’ 

2  3x  1  3x  2  3x 


+  (k.,W  -  k,.U)U  +  (k,0W  -  k. 


a  1  (h  .  3  k13U  ~  k31W)  "  h 
3x  3  3x 


3  v 
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+  a  W}  =  0 


(A- 78) 
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^  \  h!  3k1’  hl  3X1  '  \h23x2’  h23k2 


i_  1_  t,  1 _ 3) 

+  (W  “  a  h  2)  h.  *  : 

k  3  3x  3  3x 
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(A- 79) 


3  v  3  v 
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c  1  3x  1  3x  e  2  3x  2  3x 


1_  1 _ t.  1 _ 3_e  _ 1_  n2.  „  £  r 

(W  a  h  3*  h-  .  3  ”  R  el  k  G 

e  3  3x  3  3x  eff 


+  ce 2  k  ~  0 
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It  is  convenient  to  rewrite  the  governing  transport  equa¬ 
tions  ( A—  76 )  through  (A-80)  in  the  following  general  form 


V2*  =  Reff  [<a*u  -  FT  7^1  +  (' 

1  3x  1  3x 


1±  3vtk  1_  M_ 
'  h2  3X2>  h2  3X2 


+  <a*w  -  h1  — I’  S - 4  +  a*  Tt1  +  SI 

9  n3  3x  3  3x  9 


( A- 8 1 ) 


where  $  may  represent  any  one  of  the  convective  transport  quanti¬ 
ties  U,  V,  w,  k  or  e.  The  cooresponding  coefficients 
a*,  b.f  cA  and  dA  are 


=  1 f  by  2 f 

~  1 >  by  =  1» 

=  1 '  bW  =  1 ' 

=  V  bR  =  1, 
=  V  b£  =  1, 


:u  "  du  *  1 

'v  =  ^ '  dv  =  ^ 

:w  =  if  dw  »  2 

Ck  =  U  dk  =  1 
ce  =  i,  de  =  1 


( A- 8  2  ) 
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and  the  source  functions  s.  for  U,  V,  W,  k  and  e  are,  respec- 
t ively 


SU  ~  Ref f  f(k12U  "  k21V)V  +  (k13U  ~  k3iw)w  +  h1  ^1  +  3  h1  3j{l 


3  v 


h  a  1  (2k12V  +  2kl3W)  h ,  a  2  (h  1  k12U 
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3  v 
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3  3x  1  3x  2  3x 


-  2k12  hT  77  *  2k3 1  T7  7!  -  2k13  t  71  -  “llu  -  “12V 


-  a  W 
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( A-8  3a ) 
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3  v 
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-  2k23  h  7 f  +  2k  1 2  h  77  -  2k2i  h  71  -  “21“ 
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22  23 
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SW  =  Ref f  [(k31W  “  kl3LJ)CJ  +  (k32W  k23V)V  +  h3  3x3  +  3  h3  3x3 

-  r - f  (r - ^  -  k,3U  -  k3,W)  -  £ - \  - -  k23V 

hl  3X1  h3  3x 3  13  31  h2  3x2  h3  3x3  23 

3  v 

-  k32W)  -  <2k31U  +  2k32V)1  +  2k  1 3 


-  2k,,  h - Nr  +  2k,,  -r - QL  _  2k,,  r - Nf  -  a  u  -  a  V 

31  h3  3x3  23  h2  3x2  32  h3  3x3  31  32 


-  a^W 


( A-83c ) 


Sk  =  -  °k  Reff  (G  “  £) 


(A-83d) 


se  =  “  ffeReff(cel  k  G  “  Ce2  k  } 


( A-83e ) 


The  general  convective  transport  equations  (A-81)  governing 
the  transport  of  momentum  and  turbulence  quantities  can  again  be 
transformed  from  the  orthogonal  x1  coordinates  to  the  general 
curvilinear  C1  coordinates  by  using  the  expressions  for  the  gra¬ 
dient  and  Laplacian,  equations  (A-52)  and  (A-59),  respectively, 
to  obtain 


11  22  33  12  13  23  1  2 

a  ♦«+  g  *nn+  ®  2g  ♦€n+  29  »U+  29  *nC+  f  V  f  *n 


R  b 

*  cV  -r1  "v  -  /  (blvt,5k  bl\,s 


+  b2<t>n  +  b3$c)  +  [a^V - 1  (b2vt,C+  b2vt,n+  b2vt,c)J 
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(b2V  b2V  b2^}  +  [a*W  " 


-*■  (b^v  +  b,v 
J  3  t ,  4  3  t ,  n 


+  b3”t,S),(b3*E  +  b3V  b3*C>>  +  S  Ref f ^t+  s, 


( A-84 ) 


where 


SU  Reff{(k12U  "  k21V)V  +  (k13U  “  k3iw)w  +  j(biP5+  biPr,+  biPc> 
+  3J  (blk4+  blkn+  blk?)  "  j  (bivt,4+  blvt,n+  blvt,c) 
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J  2  4 


■12  ,.l 


+  b2Vn+  b2V  -  2  ~r  (blV  blV  blV  +  2  -f 
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( A-85a) 
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+  k 32V1  +  tj  (b^-t-  bxVn+  bxVc+  ^2U^+  b2Un+  b^)  -  k12U 

-  k21V]2t  >3  <bX  +  blWn  -  blV  b3b£  +  b3>V  b3Uc 1  ‘  k13U 
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It  is  noted  that  equations  (A-84)  depend  upon  g1-',  b^  , 

k  ■  •  .  a.,  and  the  Jacobian.  All  these  quantities  are  functions 

i  J  l  j 

of  coordinates  only.  When  an  analytic  transformation  is  used  to 
generate  the  coordinates,  the  geometrical  coefficients  are  known 
analytically  without  additional  approximation.  On  the  other 
hand,  if  equations  (A-60)  are  solved  numerically  for  the  coordi¬ 
nate  values  x1  =  x1  K,n,C),  then  the  coefficients  g1-1,  b^  and 
the  Jacobian  J  must  be  evaluated  numerically.  In  the  present 
study,  conventional  finite-difference  techniques  are  employed  for 
the  numerical  grid  generation  and  also  the  evaluation  of  the 
coe  fficients. 
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